From 95b5cc951d12d413331affd545c79dc49e9d159c Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Sun, 14 Dec 2025 19:48:00 -0800 Subject: [PATCH 01/22] add changes made so far --- mcdc/constant.py | 12 + mcdc/object_/element.py | 535 +++++++++++++++++++ mcdc/object_/reaction.py | 133 +++++ mcdc/transport/physics/electron/__init__.py | 3 + mcdc/transport/physics/electron/interface.py | 40 ++ mcdc/transport/physics/electron/native.py | 154 ++++++ 6 files changed, 877 insertions(+) create mode 100644 mcdc/object_/element.py create mode 100644 mcdc/transport/physics/electron/__init__.py create mode 100644 mcdc/transport/physics/electron/interface.py create mode 100644 mcdc/transport/physics/electron/native.py diff --git a/mcdc/constant.py b/mcdc/constant.py index 695ee4a5..9fa2e50e 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -27,6 +27,7 @@ SCORE_SPACE_MOMENT_CURRENT = 10 SCORE_TIME_MOMENT_MU_SQ = 11 SCORE_SPACE_MOMENT_MU_SQ = 12 +SCORE_EDEP = 13 # Tally multipliers MULTIPLIER_ENERGY = 0 @@ -96,8 +97,16 @@ REACTION_NEUTRON_FISSION_PROMPT = 5 REACTION_NEUTRON_FISSION_DELAYED = 6 +REACTION_ELECTRON_ELASTIC_SCATTERING = 101 +REACTION_ELECTRON_ELASTIC_SMALL_ANGLE = 102 +REACTION_ELECTRON_ELASTIC_LARGE_ANGLE = 103 +REACTION_ELECTRON_IONIZATION = 104 +REACTION_ELECTRON_BREMSSTRAHLUNG = 105 +REACTION_ELECTRON_EXCITATION = 106 + # Particle types PARTICLE_NEUTRON = 0 +PARTICLE_ELECTRON = 1 # Data DATA_NONE = 0 @@ -169,6 +178,9 @@ BOLTZMANN_K = 8.61733326e-5 # eV/K THERMAL_THRESHOLD_FACTOR = 400 +ELECTRON_MASS = 510.99895069e3 # eV/c^2 +ELECTRON_CUTOFF_ENERGY = 100 # eV + # Weight Windows Methods WW_USER = 0 WW_PREVIOUS = 1 diff --git a/mcdc/object_/element.py b/mcdc/object_/element.py new file mode 100644 index 00000000..adb255d8 --- /dev/null +++ b/mcdc/object_/element.py @@ -0,0 +1,535 @@ +import h5py +import numpy as np +import os + +#### + +from mcdc.reaction import ( + ReactionElectronBremsstrahlung, + ReactionElectronExcitation, + ReactionElectronElasticScattering, + ReactionElectronIonization, + decode_type, +) +from mcdc.objects import ObjectNonSingleton +from mcdc.prints import print_1d_array + +# ====================================================================================== +# Element class +# ====================================================================================== + + +class Element(ObjectNonSingleton): + def __init__(self, element_name): + label = "element" + super().__init__(label) + + # Set attributes from the hdf5 file + dir_name = os.getenv("MCDC_XSLIB") + file_name = f"{element_name}.h5" + with h5py.File(f"{dir_name}/{file_name}", "r") as f: + self.name = f["element_name"][()].decode() + self.atomic_number = f["atomic_number"][()] + self.atomic_weight_ratio = f["atomic_weight_ratio"][()] + self.xs_energy_grid = f["electron_reactions/xs_energy_grid"][()] + + self.reactions = [] + self.total_xs = np.zeros_like(self.xs_energy_grid) + + for reaction_type in f["electron_reactions"]: + if reaction_type in ["xs_energy_grid", "total"]: + continue + + if reaction_type == "bremsstrahlung": + ReactionClass = ReactionElectronBremsstrahlung + + elif reaction_type == "excitation": + ReactionClass = ReactionElectronExcitation + + elif reaction_type == "elastic_scattering": + ReactionClass = ReactionElectronElasticScattering + + elif reaction_type == "ionization": + ReactionClass = ReactionElectronIonization + + reaction = ReactionClass.from_h5_group(f[f"electron_reactions/{reaction_type}"]) + self.reactions.append(reaction) + + # Accumulate total XS + self.total_xs += reaction.xs + + def __repr__(self): + text = "\n" + text += "Element\n" + text += f" - ID: {self.ID}\n" + text += f" - Name: {self.name}\n" + text += f" - Atomic number: {self.atomic_number}\n" + text += f" - Atomic weight ratio: {self.atomic_weight_ratio}\n" + text += f" - XS energy grid {print_1d_array(self.xs_energy_grid)} eV\n" + text += " - Reactions\n" + for reaction in self.reactions: + text += f" - {decode_type(reaction.type)}\n" + return text + + +# Natural isotopic abundance data from +# https://www.nndc.bnl.gov/walletcards/search.html + +ISOTOPIC_ABUNDANCE = { + "H": { + "H1": (99.972 + 99.999) / 2, + "H2": (0.001 + 0.028) / 2, + }, + "He": { + "He3": 0.0002, + "He4": 99.9998, + }, + "Li": { + "Li6": (1.9 + 7.8) / 2, + "Li7": (92.2 + 98.1) / 2, + }, + "Be": { + "Be9": 100.0, + }, + "B": { + "B10": (18.9 + 20.4) / 2, + "B11": (79.6 + 81.1) / 2, + }, + "C": { + "C12": (98.84 + 99.04) / 2, + "C13": (0.96 + 1.16) / 2, + }, + "N": { + "N14": (99.578 + 99.663) / 2, + "N15": (0.337 + 0.422) / 2, + }, + "O": { + "O16": (99.738 + 99.776) / 2, + "O17": (0.0367 + 0.04) / 2, + "O18": (0.187 + 0.222) / 2, + }, + "F": { + "F19": 100.0, + }, + "Ne": { + "Ne20": 90.48, + "Ne21": 0.27, + "Ne22": 9.25, + }, + "Na": { + "Na23": 100.0, + }, + "Mg": { + "Mg24": (78.88 + 79.05) / 2, + "Mg25": (9.988 + 10.034) / 2, + "Mg26": (10.96 + 11.09) / 2, + }, + "Al": { + "Al27": 100.0, + }, + "Si": { + "Si28": (92.191 + 92.318) / 2, + "Si29": (4.645 + 4.699) / 2, + "Si30": (3.037 + 3.11) / 2, + }, + "P": { + "P31": 100.0, + }, + "S": { + "S32": (94.41 + 95.29) / 2, + "S33": (0.729 + 0.797) / 2, + "S34": (3.96 + 4.77) / 2, + "S36": (0.0129 + 0.0187) / 2, + }, + "Cl": { + "Cl35": (75.5 + 76.1) / 2, + "Cl37": (23.9 + 24.5) / 2, + }, + "Ar": { + "Ar36": 0.3336, + "Ar38": 0.0629, + "Ar40": 99.6035, + }, + "K": { + "K39": 93.2581, + "K40": 0.0117, + "K41": 6.7302, + }, + "Ca": { + "Ca40": 96.941, + "Ca42": 0.647, + "Ca43": 0.135, + "Ca44": 2.086, + "Ca46": 0.004, + "Ca48": 0.187, + }, + "Sc": { + "Sc45": 100.0, + }, + "Ti": { + "Ti46": 8.25, + "Ti47": 7.44, + "Ti48": 73.72, + "Ti49": 5.41, + "Ti50": 5.18, + }, + "V": { + "V50": 0.25, + "V51": 99.75, + }, + "Cr": { + "Cr50": 4.345, + "Cr52": 83.789, + "Cr53": 9.501, + "Cr54": 2.365, + }, + "Mn": { + "Mn55": 100.0, + }, + "Fe": { + "Fe54": 5.845, + "Fe56": 91.754, + "Fe57": 2.119, + "Fe58": 0.282, + }, + "Co": { + "Co59": 100.0, + }, + "Ni": { + "Ni58": 68.0769, + "Ni60": 26.2231, + "Ni61": 1.1399, + "Ni62": 3.6345, + "Ni64": 0.9256, + }, + "Cu": { + "Cu63": 69.15, + "Cu65": 30.85, + }, + "Zn": { + "Zn64": 49.17, + "Zn66": 27.73, + "Zn67": 4.04, + "Zn68": 18.45, + "Zn70": 0.61, + }, + "Ga": { + "Ga69": 60.108, + "Ga71": 39.892, + }, + "Ge": { + "Ge70": 20.52, + "Ge72": 27.45, + "Ge73": 7.76, + "Ge74": 36.52, + "Ge76": 7.75, + }, + "As": { + "As75": 100.0, + }, + "Se": { + "Se74": 0.86, + "Se76": 9.23, + "Se77": 7.6, + "Se78": 23.69, + "Se80": 49.8, + "Se82": 8.82, + }, + "Br": { + "Br79": (50.5 + 50.8) / 2, + "Br81": (49.2 + 49.5) / 2, + }, + "Kr": { + "Kr78": 0.355, + "Kr80": 2.286, + "Kr82": 11.593, + "Kr83": 11.5, + "Kr84": 56.987, + "Kr86": 17.279, + }, + "Rb": { + "Rb85": 72.17, + "Rb87": 27.83, + }, + "Sr": { + "Sr84": 0.56, + "Sr86": 9.86, + "Sr87": 7.0, + "Sr88": 82.58, + }, + "Y": { + "Y89": 100.0, + }, + "Zr": { + "Zr90": 51.45, + "Zr91": 11.22, + "Zr92": 17.15, + "Zr94": 17.38, + "Zr96": 2.8, + }, + "Nb": { + "Nb93": 100.0, + }, + "Mo": { + "Mo92": 14.649, + "Mo94": 9.187, + "Mo95": 15.873, + "Mo96": 16.673, + "Mo97": 9.582, + "Mo98": 24.292, + "Mo100": 9.744, + }, + "Ru": { + "Ru96": 5.54, + "Ru98": 1.87, + "Ru99": 12.76, + "Ru100": 12.6, + "Ru101": 17.06, + "Ru102": 31.55, + "Ru104": 18.62, + }, + "Rh": { + "Rh103": 100.0, + }, + "Pd": { + "Pd102": 1.02, + "Pd104": 11.14, + "Pd105": 22.33, + "Pd106": 27.33, + "Pd108": 26.46, + "Pd110": 11.72, + }, + "Ag": { + "Ag107": 51.839, + "Ag109": 48.161, + }, + "Cd": { + "Cd106": 1.245, + "Cd108": 0.888, + "Cd110": 12.47, + "Cd111": 12.795, + "Cd112": 24.109, + "Cd113": 12.227, + "Cd114": 28.754, + "Cd116": 7.512, + }, + "In": { + "In113": 4.281, + "In115": 95.719, + }, + "Sn": { + "Sn112": 0.97, + "Sn114": 0.66, + "Sn115": 0.34, + "Sn116": 14.54, + "Sn117": 7.68, + "Sn118": 24.22, + "Sn119": 8.59, + "Sn120": 32.58, + "Sn122": 4.63, + "Sn124": 5.79, + }, + "Sb": { + "Sb121": 57.21, + "Sb123": 42.79, + }, + "Te": { + "Te120": 0.09, + "Te122": 2.55, + "Te123": 0.89, + "Te124": 4.74, + "Te125": 7.07, + "Te126": 8.84, + "Te128": 31.74, + "Te130": 34.08, + }, + "I": { + "I127": 100.0, + }, + "Xe": { + "Xe124": 0.095, + "Xe126": 0.089, + "Xe128": 1.91, + "Xe129": 26.4, + "Xe130": 4.071, + "Xe131": 21.232, + "Xe132": 26.909, + "Xe134": 10.436, + "Xe136": 8.857, + }, + "Cs": { + "Cs133": 100.0, + }, + "Ba": { + "Ba130": 0.11, + "Ba132": 0.1, + "Ba134": 2.42, + "Ba135": 6.59, + "Ba136": 7.85, + "Ba137": 11.23, + "Ba138": 71.7, + }, + "La": { + "La138": 0.08881, + "La139": 99.91119, + }, + "Ce": { + "Ce136": 0.186, + "Ce138": 0.251, + "Ce140": 88.449, + "Ce142": 11.114, + }, + "Pr": { + "Pr141": 100.0, + }, + "Nd": { + "Nd142": 27.153, + "Nd143": 12.173, + "Nd144": 23.798, + "Nd145": 8.293, + "Nd146": 17.189, + "Nd148": 5.756, + "Nd150": 5.638, + }, + "Sm": { + "Sm144": 3.08, + "Sm147": 15.0, + "Sm149": 13.82, + "Sm150": 7.37, + "Sm152": 26.74, + "Sm154": 22.74, + }, + "Eu": { + "Eu151": 47.81, + "Eu153": 52.19, + }, + "Gd": { + "Gd152": 0.2, + "Gd154": 2.18, + "Gd155": 14.8, + "Gd156": 20.47, + "Gd157": 15.65, + "Gd158": 24.84, + "Gd160": 21.86, + }, + "Tb": { + "Tb159": 100.0, + }, + "Dy": { + "Dy156": 0.056, + "Dy158": 0.095, + "Dy160": 2.329, + "Dy161": 18.889, + "Dy162": 25.475, + "Dy163": 24.896, + "Dy164": 28.26, + }, + "Ho": { + "Ho165": 100.0, + }, + "Er": { + "Er162": 0.139, + "Er164": 1.601, + "Er166": 33.503, + "Er167": 22.869, + "Er168": 26.978, + "Er170": 14.91, + }, + "Tm": { + "Tm169": 100.0, + }, + "Yb": { + "Yb168": 0.123, + "Yb170": 2.982, + "Yb171": 14.086, + "Yb172": 21.686, + "Yb173": 16.103, + "Yb174": 32.025, + "Yb176": 12.995, + }, + "Lu": { + "Lu175": 97.401, + "Lu176": 2.599, + }, + "Hf": { + "Hf174": 0.16, + "Hf176": 5.26, + "Hf177": 18.6, + "Hf178": 27.28, + "Hf179": 13.62, + "Hf180": 35.08, + }, + "Ta": { + "Ta180": 0.01201, + "Ta181": 99.98799, + }, + "W": { + "W180": 0.12, + "W182": 26.5, + "W183": 14.31, + "W184": 30.64, + "W186": 28.43, + }, + "Re": { + "Re185": 37.4, + "Re187": 62.6, + }, + "Os": { + "Os184": 0.02, + "Os186": 1.59, + "Os187": 1.97, + "Os188": 13.24, + "Os189": 16.15, + "Os190": 26.26, + "Os192": 40.78, + }, + "Ir": { + "Ir191": 37.3, + "Ir193": 62.7, + }, + "Pt": { + "Pt190": 0.012, + "Pt192": 0.782, + "Pt194": 32.864, + "Pt195": 33.77, + "Pt196": 25.21, + "Pt198": 7.356, + }, + "Au": { + "Au197": 100.0, + }, + "Hg": { + "Hg196": 0.15, + "Hg198": 10.04, + "Hg199": 16.94, + "Hg200": 23.14, + "Hg201": 13.17, + "Hg202": 29.74, + "Hg204": 6.82, + }, + "Tl": { + "Tl203": (29.44 + 29.59) / 2, + "Tl205": (70.41 + 70.56) / 2, + }, + "Pb": { + "Pb204": 1.4, + "Pb206": 24.1, + "Pb207": 22.1, + "Pb208": 52.4, + }, + "Bi": { + "Bi209": 100.0, + }, + "Th": { + "Th230": 0.02, + "Th232": 99.98, + }, + "Pa": { + "Pa231": 100.0, + }, + "U": { + "U234": 0.0054, + "U235": 0.7204, + "U238": 99.2742, + }, +} \ No newline at end of file diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index 4b0c34c2..32cbc278 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -18,6 +18,10 @@ REACTION_NEUTRON_INELASTIC_SCATTERING, REFERENCE_FRAME_COM, REFERENCE_FRAME_LAB, + REACTION_ELECTRON_BREMSSTRAHLUNG, + REACTION_ELECTRON_EXCITATION, + REACTION_ELECTRON_ELASTIC_SCATTERING, + REACTION_ELECTRON_IONIZATION, ) from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.distribution import ( @@ -74,6 +78,14 @@ def decode_type(type_): return "Neutron inelastic scattering" elif type_ == REACTION_NEUTRON_FISSION: return "Neutron fission" + elif type_ == REACTION_ELECTRON_IONIZATION: + return "Electron ionization" + elif type_ == REACTION_ELECTRON_ELASTIC_SCATTERING: + return "Electron elastic scattering" + elif type_ == REACTION_ELECTRON_BREMSSTRAHLUNG: + return "Electron bremsstrahlung" + elif type_ == REACTION_ELECTRON_EXCITATION: + return "Electron excitation" def decode_reference_frame(type_): @@ -397,3 +409,124 @@ def set_energy_distribution(h5_group): print_error(f"Unsupported energy spectrum of type {spectrum_type}") return energy_spectrum + + +# ====================================================================================== +# Electron ionization +# ====================================================================================== + + +class ReactionElectronIonization(ReactionBase): + # Annotations for Numba mode + label: str = "electron_ionization_reaction" + # + energy_grid: NDArray[float64] + N_subshell: int + subshell_xs: list[DistributionBase] + subshell_binding_energy: NDArray[float64] + subshell_product: list[DistributionBase] + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + energy_grid, + N_subshell, + subshell_xs, + subshell_binding_energy, + subshell_product, + ): + type_ = REACTION_ELECTRON_IONIZATION + super().__init__(type_, MT, xs, xs_offset, reference_frame) + + self.energy_grid = energy_grid + self.N_subshell = N_subshell + self.subshell_xs = subshell_xs + self.subshell_binding_energy = subshell_binding_energy + self.subshell_product = subshell_product + + @classmethod + def from_h5_group(cls, h5_group): + # Basic reaction properties + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + subshells = h5_group["subshells"] + + subshell_xs = [] + subshell_binding_energy = [] + subshell_product = [] + + subshell_names = list(subshells.keys()) + + # Subshell data + for name in subshell_names: + subshell = subshells[name] + + # Subshell cross section + xs_sub = subshell["xs"][()] + energy_grid_sub = subshell["energy_grid"][()] + subshell_xs.append(DataTable(energy_grid_sub, xs_sub)) + + # Binding energy + subshell_binding_energy.append(subshell["binding_energy"][()]) + + # Secondary electron energy distribution + product_grid = subshell["product"]["energy_grid"][()] + product_offset = subshell["product"]["energy_offset"][()] + product_value = subshell["product"]["value"][()] + product_pdf = subshell["product"]["PDF"][()] + + subshell_product.append( + DataMultiPDF( + product_grid, + product_offset, + product_value, + product_pdf, + ) + ) + + N_subshell = len(subshell_names) + subshell_binding_energy = np.asarray(subshell_binding_energy) + + # Common incident energy grid (assumed identical for all subshells) + energy_grid = subshells[subshell_names[0]]["energy_grid"][()] + + return cls( + MT, + xs, + xs_offset, + reference_frame, + energy_grid, + N_subshell, + subshell_xs, + subshell_binding_energy, + subshell_product, + ) + + def __repr__(self): + text = super().__repr__() + text += f" - Electron ionization\n" + text += f" - Number of subshells: {self.N_subshell}\n" + + for i in range(self.N_subshell): + text += f" - Subshell {i + 1}\n" + text += ( + f" - Binding energy: {self.subshell_binding_energy[i]} eV\n" + ) + + xs = self.subshell_xs[i] + text += ( + f" - XS: {distribution.decode_type(xs.type)} " + f"[ID: {xs.ID}]\n" + ) + + prod = self.subshell_product[i] + text += ( + f" - Secondary electron spectrum: " + f"{distribution.decode_type(prod.type)} " + f"[ID: {prod.ID}]\n" + ) + + return text \ No newline at end of file diff --git a/mcdc/transport/physics/electron/__init__.py b/mcdc/transport/physics/electron/__init__.py new file mode 100644 index 00000000..4715cd68 --- /dev/null +++ b/mcdc/transport/physics/electron/__init__.py @@ -0,0 +1,3 @@ +import mcdc.physics.electron.interface as interface + +import mcdc.physics.electron.native as native \ No newline at end of file diff --git a/mcdc/transport/physics/electron/interface.py b/mcdc/transport/physics/electron/interface.py new file mode 100644 index 00000000..f489dc1a --- /dev/null +++ b/mcdc/transport/physics/electron/interface.py @@ -0,0 +1,40 @@ +from numba import njit + +#### + +import mcdc.transport.physics.electron.native as native + + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container, material, data): + return native.particle_speed(particle_container) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, material, particle_container, mcdc, data): + return native.macro_xs(reaction_type, material, particle_container, mcdc, data) + + +@njit +def electron_production_xs(reaction_type, particle_container, mcdc, data): + return native.electron_production_xs(reaction_type, particle_container, mcdc, data) + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, prog, data): + return native.collision(particle_container, prog, data) \ No newline at end of file diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py new file mode 100644 index 00000000..73da89d7 --- /dev/null +++ b/mcdc/transport/physics/electron/native.py @@ -0,0 +1,154 @@ +import math +import numpy as np + +from numba import njit + +#### + +import mcdc.code_factory.adapt as adapt +import mcdc.mcdc_get as mcdc_get +import mcdc.numba_types as type_ +from mcdc.print_ import print_structure +import mcdc.transport.particle as particle_module +import mcdc.transport.particle_bank as particle_bank_module +import mcdc.transport.rng as rng + +from mcdc.constant import ( + ANGLE_DISTRIBUTED, + ANGLE_ENERGY_CORRELATED, + ANGLE_ISOTROPIC, + LIGHT_SPEED, + PI, + PI_HALF, + PI_SQRT, + ELECTRON_MASS, + REACTION_TOTAL, + REACTION_ELECTRON_IONIZATION, + REACTION_ELECTRON_ELASTIC_SCATTERING, + REACTION_ELECTRON_BREMSSTRAHLUNG, + REACTION_ELECTRON_EXCITATION, +) +from mcdc.transport.data import evaluate_data +from mcdc.transport.distribution import ( + sample_correlated_distribution, + sample_distribution, + sample_isotropic_cosine, + sample_isotropic_direction, + sample_multi_table, +) +from mcdc.transport.physics.util import evaluate_xs_energy_grid, scatter_direction +from mcdc.transport.util import find_bin, linear_interpolation + + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container): + particle = particle_container[0] + E = particle["E"] + mass = ELECTRON_MASS + return LIGHT_SPEED * math.sqrt(E * (E + 2.0 * mass)) / (E + mass) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, mcdc, data): + particle = particle_container[0] + material = mcdc["native_materials"][particle["material_ID"]] + E = particle["E"] + + total = 0.0 + + for i in range(material["N_element"]): + element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) + element = mcdc["elements"][element_ID] + + element_density = mcdc_get.native_material.element_densities(i, material, data) + xs = total_micro_xs(reaction_type, E, element, data) + + total += element_density * xs + + return total + + +@njit +def total_micro_xs(reaction_type, E, element, data): + idx, E0, E1 = evaluate_xs_energy_grid(E, element, data) + if reaction_type == REACTION_TOTAL: + xs0 = mcdc_get.element.total_xs(idx, element, data) + xs1 = mcdc_get.element.total_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_IONIZATION: + xs0 = mcdc_get.element.ionization_xs(idx, element, data) + xs1 = mcdc_get.element.ionization_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_ELASTIC_SCATTERING: + xs0 = mcdc_get.element.elastic_xs(idx, element, data) + xs1 = mcdc_get.element.elastic_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_EXCITATION: + xs0 = mcdc_get.element.excitation_xs(idx, element, data) + xs1 = mcdc_get.element.excitation_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_BREMSSTRAHLUNG: + xs0 = mcdc_get.element.bremsstrahlung_xs(idx, element, data) + xs1 = mcdc_get.element.bremsstrahlung_xs(idx + 1, element, data) + return linear_interpolation(E, E0, E1, xs0, xs1) + + +@njit +def reaction_micro_xs(E, reaction_base, element, data): + idx, E0, E1 = evaluate_xs_energy_grid(E, element, data) + + # Apply offset + offset = reaction_base["xs_offset_"] + if idx < offset: + return 0.0 + else: + idx -= offset + + xs0 = mcdc_get.reaction.xs(idx, reaction_base, data) + xs1 = mcdc_get.reaction.xs(idx + 1, reaction_base, data) + return linear_interpolation(E, E0, E1, xs0, xs1) + + +@njit +def electron_production_xs(reaction_type, particle_container, mcdc, data): + particle = particle_container[0] + material_base = mcdc["materials"][particle["material_ID"]] + material = mcdc["native_materials"][material_base["child_ID"]] + + if reaction_type == REACTION_TOTAL: + ionization_type = REACTION_ELECTRON_IONIZATION + elastic_type = REACTION_ELECTRON_ELASTIC_SCATTERING + bremmstrahlung_type = REACTION_ELECTRON_BREMSSTRAHLUNG + excitation_type = REACTION_ELECTRON_EXCITATION + + ionization_xs = electron_production_xs(ionization_type, particle_container, mcdc, data) + elastic_xs = electron_production_xs(elastic_type, particle_container, mcdc, data) + bremmstrahlung_xs = electron_production_xs(bremmstrahlung_type, particle_container, mcdc, data) + excitation_xs = electron_production_xs(excitation_type, particle_container, mcdc, data) + return ionization_xs + elastic_xs + bremmstrahlung_xs + excitation_xs + + elif reaction_type == REACTION_ELECTRON_IONIZATION: + return 2 * macro_xs(reaction_type, particle_container, mcdc, data) + + elif reaction_type == REACTION_ELECTRON_ELASTIC_SCATTERING: + return macro_xs(reaction_type, particle_container, mcdc, data) + + elif reaction_type == REACTION_ELECTRON_BREMSSTRAHLUNG: + return macro_xs(reaction_type, particle_container, mcdc, data) + + elif reaction_type == REACTION_ELECTRON_EXCITATION: + return macro_xs(reaction_type, particle_container, mcdc, data) + + else: + return -1.0 + + +# ====================================================================================== +# Collision +# ====================================================================================== From a50f085f9bf2d75d15d36ce1749edeee21b98515 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:06:33 +0700 Subject: [PATCH 02/22] minor edit --- mcdc/constant.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/mcdc/constant.py b/mcdc/constant.py index 9fa2e50e..c0d7b655 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -96,7 +96,7 @@ REACTION_NEUTRON_FISSION = 4 REACTION_NEUTRON_FISSION_PROMPT = 5 REACTION_NEUTRON_FISSION_DELAYED = 6 - +# REACTION_ELECTRON_ELASTIC_SCATTERING = 101 REACTION_ELECTRON_ELASTIC_SMALL_ANGLE = 102 REACTION_ELECTRON_ELASTIC_LARGE_ANGLE = 103 @@ -172,18 +172,32 @@ AXIS_Z = 2 AXIS_T = 3 +# ====================================================================================== # Physics +# ====================================================================================== + +# Constants LIGHT_SPEED = 2.99792458e10 # cm/s -NEUTRON_MASS = 939.565413e6 # eV/c^2 BOLTZMANN_K = 8.61733326e-5 # eV/K + +# Adjustables THERMAL_THRESHOLD_FACTOR = 400 -ELECTRON_MASS = 510.99895069e3 # eV/c^2 -ELECTRON_CUTOFF_ENERGY = 100 # eV +# Neutron properties +NEUTRON_MASS = 939.565413e6 # eV/c^2 + +# Electron energy +ELECTRON_MASS = 510.99895069e3 # eV/c^2 +ELECTRON_CUTOFF_ENERGY = 100 # eV (TODO: Is this adjustable?) + +# ====================================================================================== +# Techniques +# ====================================================================================== # Weight Windows Methods WW_USER = 0 WW_PREVIOUS = 1 + # Weight Windows Modifications WW_MIN = 0 WW_WOLLABER = 1 From 2d1422b9786df2cd9ce97a54929d40a6d13405f9 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:09:53 +0700 Subject: [PATCH 03/22] edit variable types --- mcdc/object_/reaction.py | 244 +++++++++++++++++++-------------------- 1 file changed, 121 insertions(+), 123 deletions(-) diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index 32cbc278..605f28d0 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -4,6 +4,7 @@ #### +from mcdc.object_.data import DataTable import mcdc.object_.distribution as distribution from mcdc.constant import ( @@ -16,12 +17,12 @@ REACTION_NEUTRON_ELASTIC_SCATTERING, REACTION_NEUTRON_FISSION, REACTION_NEUTRON_INELASTIC_SCATTERING, - REFERENCE_FRAME_COM, - REFERENCE_FRAME_LAB, REACTION_ELECTRON_BREMSSTRAHLUNG, REACTION_ELECTRON_EXCITATION, REACTION_ELECTRON_ELASTIC_SCATTERING, REACTION_ELECTRON_IONIZATION, + REFERENCE_FRAME_COM, + REFERENCE_FRAME_LAB, ) from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.distribution import ( @@ -295,6 +296,124 @@ def __repr__(self): return text +# ====================================================================================== +# Electron ionization +# ====================================================================================== + + +class ReactionElectronIonization(ReactionBase): + # Annotations for Numba mode + label: str = "electron_ionization_reaction" + # + energy_grid: NDArray[float64] + N_subshell: int + subshell_xs: list[DataTable] + subshell_binding_energy: NDArray[float64] + subshell_product: list[DistributionMultiTable] + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + energy_grid, + N_subshell, + subshell_xs, + subshell_binding_energy, + subshell_product, + ): + type_ = REACTION_ELECTRON_IONIZATION + super().__init__(type_, MT, xs, xs_offset, reference_frame) + + self.energy_grid = energy_grid + self.N_subshell = N_subshell + self.subshell_xs = subshell_xs + self.subshell_binding_energy = subshell_binding_energy + self.subshell_product = subshell_product + + @classmethod + def from_h5_group(cls, h5_group): + # Basic reaction properties + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + subshells = h5_group["subshells"] + + subshell_xs = [] + subshell_binding_energy = [] + subshell_product = [] + + subshell_names = list(subshells.keys()) + + # Subshell data + for name in subshell_names: + subshell = subshells[name] + + # Subshell cross section + xs_sub = subshell["xs"][()] + energy_grid_sub = subshell["energy_grid"][()] + subshell_xs.append(DataTable(energy_grid_sub, xs_sub)) + + # Binding energy + subshell_binding_energy.append(subshell["binding_energy"][()]) + + # Secondary electron energy distribution + product_grid = subshell["product"]["energy_grid"][()] + product_offset = subshell["product"]["energy_offset"][()] + product_value = subshell["product"]["value"][()] + product_pdf = subshell["product"]["PDF"][()] + + subshell_product.append( + DistributionMultiTable( + product_grid, + product_offset, + product_value, + product_pdf, + ) + ) + + N_subshell = len(subshell_names) + subshell_binding_energy = subshell_binding_energy + + # Common incident energy grid (assumed identical for all subshells) + energy_grid = subshells[subshell_names[0]]["energy_grid"][()] + + return cls( + MT, + xs, + xs_offset, + reference_frame, + energy_grid, + N_subshell, + subshell_xs, + subshell_binding_energy, + subshell_product, + ) + + def __repr__(self): + text = super().__repr__() + text += f" - Electron ionization\n" + text += f" - Number of subshells: {self.N_subshell}\n" + + for i in range(self.N_subshell): + text += f" - Subshell {i + 1}\n" + text += f" - Binding energy: {self.subshell_binding_energy[i]} eV\n" + + xs = self.subshell_xs[i] + text += ( + f" - XS: {distribution.decode_type(xs.type)} " f"[ID: {xs.ID}]\n" + ) + + prod = self.subshell_product[i] + text += ( + f" - Secondary electron spectrum: " + f"{distribution.decode_type(prod.type)} " + f"[ID: {prod.ID}]\n" + ) + + return text + + # ====================================================================================== # Helper functions # ====================================================================================== @@ -409,124 +528,3 @@ def set_energy_distribution(h5_group): print_error(f"Unsupported energy spectrum of type {spectrum_type}") return energy_spectrum - - -# ====================================================================================== -# Electron ionization -# ====================================================================================== - - -class ReactionElectronIonization(ReactionBase): - # Annotations for Numba mode - label: str = "electron_ionization_reaction" - # - energy_grid: NDArray[float64] - N_subshell: int - subshell_xs: list[DistributionBase] - subshell_binding_energy: NDArray[float64] - subshell_product: list[DistributionBase] - - def __init__( - self, - MT, - xs, - xs_offset, - reference_frame, - energy_grid, - N_subshell, - subshell_xs, - subshell_binding_energy, - subshell_product, - ): - type_ = REACTION_ELECTRON_IONIZATION - super().__init__(type_, MT, xs, xs_offset, reference_frame) - - self.energy_grid = energy_grid - self.N_subshell = N_subshell - self.subshell_xs = subshell_xs - self.subshell_binding_energy = subshell_binding_energy - self.subshell_product = subshell_product - - @classmethod - def from_h5_group(cls, h5_group): - # Basic reaction properties - MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) - - subshells = h5_group["subshells"] - - subshell_xs = [] - subshell_binding_energy = [] - subshell_product = [] - - subshell_names = list(subshells.keys()) - - # Subshell data - for name in subshell_names: - subshell = subshells[name] - - # Subshell cross section - xs_sub = subshell["xs"][()] - energy_grid_sub = subshell["energy_grid"][()] - subshell_xs.append(DataTable(energy_grid_sub, xs_sub)) - - # Binding energy - subshell_binding_energy.append(subshell["binding_energy"][()]) - - # Secondary electron energy distribution - product_grid = subshell["product"]["energy_grid"][()] - product_offset = subshell["product"]["energy_offset"][()] - product_value = subshell["product"]["value"][()] - product_pdf = subshell["product"]["PDF"][()] - - subshell_product.append( - DataMultiPDF( - product_grid, - product_offset, - product_value, - product_pdf, - ) - ) - - N_subshell = len(subshell_names) - subshell_binding_energy = np.asarray(subshell_binding_energy) - - # Common incident energy grid (assumed identical for all subshells) - energy_grid = subshells[subshell_names[0]]["energy_grid"][()] - - return cls( - MT, - xs, - xs_offset, - reference_frame, - energy_grid, - N_subshell, - subshell_xs, - subshell_binding_energy, - subshell_product, - ) - - def __repr__(self): - text = super().__repr__() - text += f" - Electron ionization\n" - text += f" - Number of subshells: {self.N_subshell}\n" - - for i in range(self.N_subshell): - text += f" - Subshell {i + 1}\n" - text += ( - f" - Binding energy: {self.subshell_binding_energy[i]} eV\n" - ) - - xs = self.subshell_xs[i] - text += ( - f" - XS: {distribution.decode_type(xs.type)} " - f"[ID: {xs.ID}]\n" - ) - - prod = self.subshell_product[i] - text += ( - f" - Secondary electron spectrum: " - f"{distribution.decode_type(prod.type)} " - f"[ID: {prod.ID}]\n" - ) - - return text \ No newline at end of file From 9a7707c6ef9d5eee25f4a820b0d7f86370d852bd Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:17:13 +0700 Subject: [PATCH 04/22] add elements to simulation object --- mcdc/object_/simulation.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py index 89076d86..6383e307 100644 --- a/mcdc/object_/simulation.py +++ b/mcdc/object_/simulation.py @@ -10,6 +10,7 @@ if TYPE_CHECKING: from mcdc.object_.cell import Cell, Region + from mcdc.object_.element import Element from mcdc.object_.material import MaterialBase from mcdc.object_.nuclide import Nuclide from mcdc.object_.reaction import ReactionBase @@ -57,6 +58,7 @@ class Simulation(ObjectSingleton): data: list[DataBase] distributions: list[DistributionBase] materials: list[MaterialBase] + elements: list[Element] nuclides: list[Nuclide] reactions: list[ReactionBase] sources: list[Source] @@ -139,6 +141,7 @@ def __init__(self): self.data = [DataNone()] self.distributions = [DistributionNone()] self.materials = [] + self.elements = [] self.nuclides = [] self.reactions = [] self.sources = [] From d9509a2afd52a831963a82f304e2a80a4ce6cd19 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:17:42 +0700 Subject: [PATCH 05/22] minor edit --- mcdc/object_/element.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/mcdc/object_/element.py b/mcdc/object_/element.py index adb255d8..a0a60da5 100644 --- a/mcdc/object_/element.py +++ b/mcdc/object_/element.py @@ -4,15 +4,15 @@ #### -from mcdc.reaction import ( +from mcdc.object_.reaction import ( ReactionElectronBremsstrahlung, ReactionElectronExcitation, ReactionElectronElasticScattering, ReactionElectronIonization, decode_type, ) -from mcdc.objects import ObjectNonSingleton -from mcdc.prints import print_1d_array +from mcdc.object_.base import ObjectNonSingleton +from mcdc.print_ import print_1d_array # ====================================================================================== # Element class @@ -20,6 +20,10 @@ class Element(ObjectNonSingleton): + # Annotations for Numba mode + label: str = "element" + # + def __init__(self, element_name): label = "element" super().__init__(label) @@ -52,7 +56,9 @@ def __init__(self, element_name): elif reaction_type == "ionization": ReactionClass = ReactionElectronIonization - reaction = ReactionClass.from_h5_group(f[f"electron_reactions/{reaction_type}"]) + reaction = ReactionClass.from_h5_group( + f[f"electron_reactions/{reaction_type}"] + ) self.reactions.append(reaction) # Accumulate total XS @@ -72,7 +78,7 @@ def __repr__(self): return text -# Natural isotopic abundance data from +# Natural isotopic abundance data from # https://www.nndc.bnl.gov/walletcards/search.html ISOTOPIC_ABUNDANCE = { @@ -532,4 +538,4 @@ def __repr__(self): "U235": 0.7204, "U238": 99.2742, }, -} \ No newline at end of file +} From 3d336dcf6e8d52c870fa9b030eebd4108553269b Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:19:45 +0700 Subject: [PATCH 06/22] add a data placeholder in object util --- mcdc/object_/util.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/mcdc/object_/util.py b/mcdc/object_/util.py index efa3e74c..dd63685b 100644 --- a/mcdc/object_/util.py +++ b/mcdc/object_/util.py @@ -59,3 +59,8 @@ def multi_cdf_from_pdf(offset, value, pdf): def is_sorted(a): return np.all(a[:-1] <= a[1:]) + + +# ====================================================================================== +# Data +# ====================================================================================== From b40c0495312178d8fd9e1072b76a0eb4a5aef26c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Tue, 16 Dec 2025 15:28:00 +0700 Subject: [PATCH 07/22] revert back to using base type --- mcdc/object_/reaction.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index 605f28d0..b20bf68d 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -4,7 +4,7 @@ #### -from mcdc.object_.data import DataTable +from mcdc.object_.data import DataBase, DataTable import mcdc.object_.distribution as distribution from mcdc.constant import ( @@ -307,9 +307,9 @@ class ReactionElectronIonization(ReactionBase): # energy_grid: NDArray[float64] N_subshell: int - subshell_xs: list[DataTable] + subshell_xs: list[DataBase] subshell_binding_energy: NDArray[float64] - subshell_product: list[DistributionMultiTable] + subshell_product: list[DistributionBase] def __init__( self, From e00a0c3b300a592e1e64207acb141b08a5caa470 Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Wed, 17 Dec 2025 23:53:30 -0800 Subject: [PATCH 08/22] fix: correct fundamental electron transport physics to match new refactoring --- mcdc/constant.py | 2 + mcdc/object_/element.py | 661 +++++----------------- mcdc/object_/reaction.py | 203 +++++-- mcdc/object_/util.py | 462 +++++++++++++++ mcdc/transport/physics/electron/native.py | 472 ++++++++++++++- 5 files changed, 1216 insertions(+), 584 deletions(-) diff --git a/mcdc/constant.py b/mcdc/constant.py index c0d7b655..9c19ef2c 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -189,6 +189,8 @@ # Electron energy ELECTRON_MASS = 510.99895069e3 # eV/c^2 ELECTRON_CUTOFF_ENERGY = 100 # eV (TODO: Is this adjustable?) +FINE_STRUCTURE_CONSTANT = 7.2973525693e-3 # ~1/137 +MU_CUTOFF = 0.999999 # ====================================================================================== # Techniques diff --git a/mcdc/object_/element.py b/mcdc/object_/element.py index a0a60da5..5e808aab 100644 --- a/mcdc/object_/element.py +++ b/mcdc/object_/element.py @@ -2,20 +2,23 @@ import numpy as np import os +from numpy import float64 +from numpy.typing import NDArray + #### +from mcdc.object_.base import ObjectNonSingleton from mcdc.object_.reaction import ( - ReactionElectronBremsstrahlung, - ReactionElectronExcitation, ReactionElectronElasticScattering, + ReactionElectronExcitation, + ReactionElectronBremsstrahlung, ReactionElectronIonization, - decode_type, ) -from mcdc.object_.base import ObjectNonSingleton -from mcdc.print_ import print_1d_array +from mcdc.print_ import print_1d_array, print_error + # ====================================================================================== -# Element class +# Element # ====================================================================================== @@ -23,519 +26,153 @@ class Element(ObjectNonSingleton): # Annotations for Numba mode label: str = "element" # - - def __init__(self, element_name): - label = "element" - super().__init__(label) + name: str + atomic_weight_ratio: float + atomic_number: int + # XS on xs_energy_grid + xs_energy_grid: NDArray[float64] + total_xs: NDArray[float64] + ionization_xs: NDArray[float64] + elastic_xs: NDArray[float64] + excitation_xs: NDArray[float64] + bremsstrahlung_xs: NDArray[float64] + # Reactions + ionization_reactions: list[ReactionElectronIonization] + elastic_scattering_reactions: list[ReactionElectronElasticScattering] + excitation_reactions: list[ReactionElectronExcitation] + bremsstrahlung_reactions: list[ReactionElectronBremsstrahlung] + # Ionization element attributes + ionization_subshell_binding_energy: NDArray[float64] + + def __init__(self, element_name: str): + super().__init__() + + self.name = element_name # Set attributes from the hdf5 file - dir_name = os.getenv("MCDC_XSLIB") + dir_name = os.getenv("MCDC_LIB") file_name = f"{element_name}.h5" - with h5py.File(f"{dir_name}/{file_name}", "r") as f: - self.name = f["element_name"][()].decode() - self.atomic_number = f["atomic_number"][()] - self.atomic_weight_ratio = f["atomic_weight_ratio"][()] - self.xs_energy_grid = f["electron_reactions/xs_energy_grid"][()] - - self.reactions = [] - self.total_xs = np.zeros_like(self.xs_energy_grid) - - for reaction_type in f["electron_reactions"]: - if reaction_type in ["xs_energy_grid", "total"]: - continue - - if reaction_type == "bremsstrahlung": - ReactionClass = ReactionElectronBremsstrahlung - - elif reaction_type == "excitation": - ReactionClass = ReactionElectronExcitation - - elif reaction_type == "elastic_scattering": - ReactionClass = ReactionElectronElasticScattering - - elif reaction_type == "ionization": - ReactionClass = ReactionElectronIonization - - reaction = ReactionClass.from_h5_group( - f[f"electron_reactions/{reaction_type}"] - ) - self.reactions.append(reaction) - - # Accumulate total XS - self.total_xs += reaction.xs + file = h5py.File(f"{dir_name}/{file_name}", "r") + + # Basic properties + self.atomic_weight_ratio = float(file["atomic_weight_ratio"][()]) + self.atomic_number = int(file["atomic_number"][()]) + + # The reactions + rx_names = [ + "elastic_scattering", + "excitation", + "bremsstrahlung", + "ionization", + ] + + # The reaction MTs + MTs = {} + for name in rx_names: + if name not in file["electron_reactions"]: + MTs[name] = [] + continue + + MTs[name] = [ + x for x in file[f"electron_reactions/{name}"] if x.startswith("MT") + ] + + # ========================================================================== + # Reaction XS + # ========================================================================== + + # Energy grid + xs_energy = file["electron_reactions/xs_energy_grid"][()] + self.xs_energy_grid = xs_energy + + # The total XS + self.total_xs = np.zeros_like(self.xs_energy_grid) + self.elastic_xs = np.zeros_like(self.xs_energy_grid) + self.excitation_xs = np.zeros_like(self.xs_energy_grid) + self.bremsstrahlung_xs = np.zeros_like(self.xs_energy_grid) + self.ionization_xs = np.zeros_like(self.xs_energy_grid) + + xs_containers = [ + self.elastic_xs, + self.excitation_xs, + self.bremsstrahlung_xs, + self.ionization_xs, + ] + for xs_container, rx_name in list(zip(xs_containers, rx_names)): + for MT in MTs[rx_name]: + xs = file[f"electron_reactions/{rx_name}/{MT}/xs"] + xs_container[xs.attrs["offset"] :] += xs[()] + + # Total = sum of xs' + self.total_xs = ( + self.elastic_xs + + self.excitation_xs + + self.bremsstrahlung_xs + + self.ionization_xs + ) + + # ========================================================================== + # The reactions + # ========================================================================== + + self.elastic_scattering_reactions = [] + self.excitation_reactions = [] + self.bremsstrahlung_reactions = [] + self.ionization_reactions = [] + + rx_containers = [ + self.elastic_scattering_reactions, + self.excitation_reactions, + self.bremsstrahlung_reactions, + self.ionization_reactions, + ] + rx_classes = [ + ReactionElectronElasticScattering, + ReactionElectronExcitation, + ReactionElectronBremsstrahlung, + ReactionElectronIonization, + ] + + for rx_container, rx_name, rx_class in list(zip(rx_containers, rx_names, rx_classes)): + for MT in MTs[rx_name]: + h5_group = file[f"electron_reactions/{rx_name}/{MT}"] + reaction = rx_class.from_h5_group(h5_group) + rx_container.append(reaction) + + # ========================================================================== + # Ionization element attributes + # ========================================================================== + + be_list = [] + if len(MTs["ionization"]) > 0: + g = file[f"electron_reactions/ionization/{MTs['ionization'][0]}"] + subshells = g["subshells"] + subshell_names = list(subshells.keys()) + for name in subshell_names: + sh = subshells[name] + be_list.append(float(sh["binding_energy"][()])) + + self.ionization_subshell_binding_energy = np.asarray(be_list) + + file.close() def __repr__(self): text = "\n" - text += "Element\n" + text += f"Element\n" text += f" - ID: {self.ID}\n" text += f" - Name: {self.name}\n" text += f" - Atomic number: {self.atomic_number}\n" text += f" - Atomic weight ratio: {self.atomic_weight_ratio}\n" - text += f" - XS energy grid {print_1d_array(self.xs_energy_grid)} eV\n" - text += " - Reactions\n" - for reaction in self.reactions: - text += f" - {decode_type(reaction.type)}\n" - return text - - -# Natural isotopic abundance data from -# https://www.nndc.bnl.gov/walletcards/search.html - -ISOTOPIC_ABUNDANCE = { - "H": { - "H1": (99.972 + 99.999) / 2, - "H2": (0.001 + 0.028) / 2, - }, - "He": { - "He3": 0.0002, - "He4": 99.9998, - }, - "Li": { - "Li6": (1.9 + 7.8) / 2, - "Li7": (92.2 + 98.1) / 2, - }, - "Be": { - "Be9": 100.0, - }, - "B": { - "B10": (18.9 + 20.4) / 2, - "B11": (79.6 + 81.1) / 2, - }, - "C": { - "C12": (98.84 + 99.04) / 2, - "C13": (0.96 + 1.16) / 2, - }, - "N": { - "N14": (99.578 + 99.663) / 2, - "N15": (0.337 + 0.422) / 2, - }, - "O": { - "O16": (99.738 + 99.776) / 2, - "O17": (0.0367 + 0.04) / 2, - "O18": (0.187 + 0.222) / 2, - }, - "F": { - "F19": 100.0, - }, - "Ne": { - "Ne20": 90.48, - "Ne21": 0.27, - "Ne22": 9.25, - }, - "Na": { - "Na23": 100.0, - }, - "Mg": { - "Mg24": (78.88 + 79.05) / 2, - "Mg25": (9.988 + 10.034) / 2, - "Mg26": (10.96 + 11.09) / 2, - }, - "Al": { - "Al27": 100.0, - }, - "Si": { - "Si28": (92.191 + 92.318) / 2, - "Si29": (4.645 + 4.699) / 2, - "Si30": (3.037 + 3.11) / 2, - }, - "P": { - "P31": 100.0, - }, - "S": { - "S32": (94.41 + 95.29) / 2, - "S33": (0.729 + 0.797) / 2, - "S34": (3.96 + 4.77) / 2, - "S36": (0.0129 + 0.0187) / 2, - }, - "Cl": { - "Cl35": (75.5 + 76.1) / 2, - "Cl37": (23.9 + 24.5) / 2, - }, - "Ar": { - "Ar36": 0.3336, - "Ar38": 0.0629, - "Ar40": 99.6035, - }, - "K": { - "K39": 93.2581, - "K40": 0.0117, - "K41": 6.7302, - }, - "Ca": { - "Ca40": 96.941, - "Ca42": 0.647, - "Ca43": 0.135, - "Ca44": 2.086, - "Ca46": 0.004, - "Ca48": 0.187, - }, - "Sc": { - "Sc45": 100.0, - }, - "Ti": { - "Ti46": 8.25, - "Ti47": 7.44, - "Ti48": 73.72, - "Ti49": 5.41, - "Ti50": 5.18, - }, - "V": { - "V50": 0.25, - "V51": 99.75, - }, - "Cr": { - "Cr50": 4.345, - "Cr52": 83.789, - "Cr53": 9.501, - "Cr54": 2.365, - }, - "Mn": { - "Mn55": 100.0, - }, - "Fe": { - "Fe54": 5.845, - "Fe56": 91.754, - "Fe57": 2.119, - "Fe58": 0.282, - }, - "Co": { - "Co59": 100.0, - }, - "Ni": { - "Ni58": 68.0769, - "Ni60": 26.2231, - "Ni61": 1.1399, - "Ni62": 3.6345, - "Ni64": 0.9256, - }, - "Cu": { - "Cu63": 69.15, - "Cu65": 30.85, - }, - "Zn": { - "Zn64": 49.17, - "Zn66": 27.73, - "Zn67": 4.04, - "Zn68": 18.45, - "Zn70": 0.61, - }, - "Ga": { - "Ga69": 60.108, - "Ga71": 39.892, - }, - "Ge": { - "Ge70": 20.52, - "Ge72": 27.45, - "Ge73": 7.76, - "Ge74": 36.52, - "Ge76": 7.75, - }, - "As": { - "As75": 100.0, - }, - "Se": { - "Se74": 0.86, - "Se76": 9.23, - "Se77": 7.6, - "Se78": 23.69, - "Se80": 49.8, - "Se82": 8.82, - }, - "Br": { - "Br79": (50.5 + 50.8) / 2, - "Br81": (49.2 + 49.5) / 2, - }, - "Kr": { - "Kr78": 0.355, - "Kr80": 2.286, - "Kr82": 11.593, - "Kr83": 11.5, - "Kr84": 56.987, - "Kr86": 17.279, - }, - "Rb": { - "Rb85": 72.17, - "Rb87": 27.83, - }, - "Sr": { - "Sr84": 0.56, - "Sr86": 9.86, - "Sr87": 7.0, - "Sr88": 82.58, - }, - "Y": { - "Y89": 100.0, - }, - "Zr": { - "Zr90": 51.45, - "Zr91": 11.22, - "Zr92": 17.15, - "Zr94": 17.38, - "Zr96": 2.8, - }, - "Nb": { - "Nb93": 100.0, - }, - "Mo": { - "Mo92": 14.649, - "Mo94": 9.187, - "Mo95": 15.873, - "Mo96": 16.673, - "Mo97": 9.582, - "Mo98": 24.292, - "Mo100": 9.744, - }, - "Ru": { - "Ru96": 5.54, - "Ru98": 1.87, - "Ru99": 12.76, - "Ru100": 12.6, - "Ru101": 17.06, - "Ru102": 31.55, - "Ru104": 18.62, - }, - "Rh": { - "Rh103": 100.0, - }, - "Pd": { - "Pd102": 1.02, - "Pd104": 11.14, - "Pd105": 22.33, - "Pd106": 27.33, - "Pd108": 26.46, - "Pd110": 11.72, - }, - "Ag": { - "Ag107": 51.839, - "Ag109": 48.161, - }, - "Cd": { - "Cd106": 1.245, - "Cd108": 0.888, - "Cd110": 12.47, - "Cd111": 12.795, - "Cd112": 24.109, - "Cd113": 12.227, - "Cd114": 28.754, - "Cd116": 7.512, - }, - "In": { - "In113": 4.281, - "In115": 95.719, - }, - "Sn": { - "Sn112": 0.97, - "Sn114": 0.66, - "Sn115": 0.34, - "Sn116": 14.54, - "Sn117": 7.68, - "Sn118": 24.22, - "Sn119": 8.59, - "Sn120": 32.58, - "Sn122": 4.63, - "Sn124": 5.79, - }, - "Sb": { - "Sb121": 57.21, - "Sb123": 42.79, - }, - "Te": { - "Te120": 0.09, - "Te122": 2.55, - "Te123": 0.89, - "Te124": 4.74, - "Te125": 7.07, - "Te126": 8.84, - "Te128": 31.74, - "Te130": 34.08, - }, - "I": { - "I127": 100.0, - }, - "Xe": { - "Xe124": 0.095, - "Xe126": 0.089, - "Xe128": 1.91, - "Xe129": 26.4, - "Xe130": 4.071, - "Xe131": 21.232, - "Xe132": 26.909, - "Xe134": 10.436, - "Xe136": 8.857, - }, - "Cs": { - "Cs133": 100.0, - }, - "Ba": { - "Ba130": 0.11, - "Ba132": 0.1, - "Ba134": 2.42, - "Ba135": 6.59, - "Ba136": 7.85, - "Ba137": 11.23, - "Ba138": 71.7, - }, - "La": { - "La138": 0.08881, - "La139": 99.91119, - }, - "Ce": { - "Ce136": 0.186, - "Ce138": 0.251, - "Ce140": 88.449, - "Ce142": 11.114, - }, - "Pr": { - "Pr141": 100.0, - }, - "Nd": { - "Nd142": 27.153, - "Nd143": 12.173, - "Nd144": 23.798, - "Nd145": 8.293, - "Nd146": 17.189, - "Nd148": 5.756, - "Nd150": 5.638, - }, - "Sm": { - "Sm144": 3.08, - "Sm147": 15.0, - "Sm149": 13.82, - "Sm150": 7.37, - "Sm152": 26.74, - "Sm154": 22.74, - }, - "Eu": { - "Eu151": 47.81, - "Eu153": 52.19, - }, - "Gd": { - "Gd152": 0.2, - "Gd154": 2.18, - "Gd155": 14.8, - "Gd156": 20.47, - "Gd157": 15.65, - "Gd158": 24.84, - "Gd160": 21.86, - }, - "Tb": { - "Tb159": 100.0, - }, - "Dy": { - "Dy156": 0.056, - "Dy158": 0.095, - "Dy160": 2.329, - "Dy161": 18.889, - "Dy162": 25.475, - "Dy163": 24.896, - "Dy164": 28.26, - }, - "Ho": { - "Ho165": 100.0, - }, - "Er": { - "Er162": 0.139, - "Er164": 1.601, - "Er166": 33.503, - "Er167": 22.869, - "Er168": 26.978, - "Er170": 14.91, - }, - "Tm": { - "Tm169": 100.0, - }, - "Yb": { - "Yb168": 0.123, - "Yb170": 2.982, - "Yb171": 14.086, - "Yb172": 21.686, - "Yb173": 16.103, - "Yb174": 32.025, - "Yb176": 12.995, - }, - "Lu": { - "Lu175": 97.401, - "Lu176": 2.599, - }, - "Hf": { - "Hf174": 0.16, - "Hf176": 5.26, - "Hf177": 18.6, - "Hf178": 27.28, - "Hf179": 13.62, - "Hf180": 35.08, - }, - "Ta": { - "Ta180": 0.01201, - "Ta181": 99.98799, - }, - "W": { - "W180": 0.12, - "W182": 26.5, - "W183": 14.31, - "W184": 30.64, - "W186": 28.43, - }, - "Re": { - "Re185": 37.4, - "Re187": 62.6, - }, - "Os": { - "Os184": 0.02, - "Os186": 1.59, - "Os187": 1.97, - "Os188": 13.24, - "Os189": 16.15, - "Os190": 26.26, - "Os192": 40.78, - }, - "Ir": { - "Ir191": 37.3, - "Ir193": 62.7, - }, - "Pt": { - "Pt190": 0.012, - "Pt192": 0.782, - "Pt194": 32.864, - "Pt195": 33.77, - "Pt196": 25.21, - "Pt198": 7.356, - }, - "Au": { - "Au197": 100.0, - }, - "Hg": { - "Hg196": 0.15, - "Hg198": 10.04, - "Hg199": 16.94, - "Hg200": 23.14, - "Hg201": 13.17, - "Hg202": 29.74, - "Hg204": 6.82, - }, - "Tl": { - "Tl203": (29.44 + 29.59) / 2, - "Tl205": (70.41 + 70.56) / 2, - }, - "Pb": { - "Pb204": 1.4, - "Pb206": 24.1, - "Pb207": 22.1, - "Pb208": 52.4, - }, - "Bi": { - "Bi209": 100.0, - }, - "Th": { - "Th230": 0.02, - "Th232": 99.98, - }, - "Pa": { - "Pa231": 100.0, - }, - "U": { - "U234": 0.0054, - "U235": 0.7204, - "U238": 99.2742, - }, -} + text += f" - Reaction MTs\n" + text += f" - Ionization: {[int(x.MT) for x in self.ionization_reactions]}\n" + text += f" - Elastic scattering: {[int(x.MT) for x in self.elastic_scattering_reactions]}\n" + text += f" - Excitation: {[int(x.MT) for x in self.excitation_reactions]}\n" + text += f" - Bremsstrahlung: {[int(x.MT) for x in self.bremsstrahlung_reactions]}\n" + text += f" - Reaction cross-sections (eV, barns)\n" + text += f" - Energy grid {print_1d_array(self.xs_energy_grid)}\n" + text += f" - Total {print_1d_array(self.total_xs)}\n" + text += f" - Ionization {print_1d_array(self.ionization_xs)}\n" + text += f" - Elastic scattering {print_1d_array(self.elastic_xs)}\n" + text += f" - Excitation {print_1d_array(self.excitation_xs)}\n" + text += f" - Bremsstrahlung {print_1d_array(self.bremsstrahlung_xs)}\n" + return text \ No newline at end of file diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index b20bf68d..2a9684c2 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -302,13 +302,10 @@ def __repr__(self): class ReactionElectronIonization(ReactionBase): - # Annotations for Numba mode label: str = "electron_ionization_reaction" # - energy_grid: NDArray[float64] N_subshell: int subshell_xs: list[DataBase] - subshell_binding_energy: NDArray[float64] subshell_product: list[DistributionBase] def __init__( @@ -317,100 +314,196 @@ def __init__( xs, xs_offset, reference_frame, - energy_grid, - N_subshell, subshell_xs, - subshell_binding_energy, subshell_product, ): type_ = REACTION_ELECTRON_IONIZATION super().__init__(type_, MT, xs, xs_offset, reference_frame) - self.energy_grid = energy_grid - self.N_subshell = N_subshell + self.N_subshell = len(subshell_xs) self.subshell_xs = subshell_xs - self.subshell_binding_energy = subshell_binding_energy self.subshell_product = subshell_product @classmethod def from_h5_group(cls, h5_group): - # Basic reaction properties MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) subshells = h5_group["subshells"] + subshell_names = list(subshells.keys()) subshell_xs = [] - subshell_binding_energy = [] subshell_product = [] - subshell_names = list(subshells.keys()) - - # Subshell data for name in subshell_names: subshell = subshells[name] - # Subshell cross section - xs_sub = subshell["xs"][()] - energy_grid_sub = subshell["energy_grid"][()] - subshell_xs.append(DataTable(energy_grid_sub, xs_sub)) - - # Binding energy - subshell_binding_energy.append(subshell["binding_energy"][()]) + # Subshell cross section table (each has its own energy grid) + xs_sub = DataTable( + subshell["energy_grid"][()], + subshell["xs"][()] + ) + subshell_xs.append(xs_sub) # Secondary electron energy distribution - product_grid = subshell["product"]["energy_grid"][()] - product_offset = subshell["product"]["energy_offset"][()] - product_value = subshell["product"]["value"][()] - product_pdf = subshell["product"]["PDF"][()] - - subshell_product.append( - DistributionMultiTable( - product_grid, - product_offset, - product_value, - product_pdf, - ) + product = subshell["product"] + prod_dist = DistributionMultiTable( + product["energy_grid"][()], + product["energy_offset"][()], + product["value"][()], + product["PDF"][()] ) - - N_subshell = len(subshell_names) - subshell_binding_energy = subshell_binding_energy - - # Common incident energy grid (assumed identical for all subshells) - energy_grid = subshells[subshell_names[0]]["energy_grid"][()] + subshell_product.append(prod_dist) return cls( MT, xs, xs_offset, reference_frame, - energy_grid, - N_subshell, subshell_xs, - subshell_binding_energy, subshell_product, ) def __repr__(self): text = super().__repr__() - text += f" - Electron ionization\n" text += f" - Number of subshells: {self.N_subshell}\n" - for i in range(self.N_subshell): - text += f" - Subshell {i + 1}\n" - text += f" - Binding energy: {self.subshell_binding_energy[i]} eV\n" - - xs = self.subshell_xs[i] - text += ( - f" - XS: {distribution.decode_type(xs.type)} " f"[ID: {xs.ID}]\n" - ) - + text += f" - Subshell {i+1}\n" + text += f" - XS: DataTable [ID: {self.subshell_xs[i].ID}]\n" prod = self.subshell_product[i] text += ( f" - Secondary electron spectrum: " - f"{distribution.decode_type(prod.type)} " - f"[ID: {prod.ID}]\n" + f"{distribution.decode_type(prod.type)} [ID: {prod.ID}]\n" ) + return text + + +# ====================================================================================== +# Electron elastic scattering +# ====================================================================================== + + +class ReactionElectronElasticScattering(ReactionBase): + label: str = "electron_elastic_scattering_reaction" + # + mu_cut: float + xs_large: DataBase + mu_table: DistributionMultiTable + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + mu_cut, + xs_large, + mu_table, + ): + type_ = REACTION_ELECTRON_ELASTIC_SCATTERING + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.mu_cut = mu_cut + self.xs_large = xs_large + self.mu_table = mu_table + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + # Large angle data + large_angle = h5_group["large_angle"] + mu_cut = float(large_angle.attrs.get("mu_cut", 0.999999)) + + xs_large = DataTable( + large_angle["xs_energy"][()], + large_angle["xs"][()] + ) + + mu_g = large_angle["scattering_cosine"] + mu_table = DistributionMultiTable( + mu_g["energy_grid"][()], + mu_g["energy_offset"][()], + mu_g["value"][()], + mu_g["PDF"][()] + ) + + return cls( + MT, + xs, + xs_offset, + reference_frame, + mu_cut, + xs_large, + mu_table, + ) + + def __repr__(self): + text = super().__repr__() + text += f" - Mu cut: {self.mu_cut}\n" + text += f" - Large angle XS: DataTable [ID: {self.xs_large.ID}]\n" + text += ( + f" - Scattering cosine: {distribution.decode_type(self.mu_table.type)} " + f"[ID: {self.mu_table.ID}]\n" + ) + return text + + +# ====================================================================================== +# Electron bremsstrahlung +# ====================================================================================== + + +class ReactionElectronBremsstrahlung(ReactionBase): + label: str = "electron_bremsstrahlung_reaction" + # + eloss: DataBase + + def __init__(self, MT, xs, xs_offset, reference_frame, eloss): + type_ = REACTION_ELECTRON_BREMSSTRAHLUNG + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.eloss = eloss + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + base = h5_group["energy_loss"] + eloss = DataTable(base["energy"][()], base["value"][()]) + return cls(MT, xs, xs_offset, reference_frame, eloss) + + def __repr__(self): + text = super().__repr__() + text += f" - Energy loss: DataTable [ID: {self.eloss.ID}]\n" + return text + + +# ====================================================================================== +# Electron excitation +# ====================================================================================== + + +class ReactionElectronExcitation(ReactionBase): + label: str = "electron_excitation_reaction" + # + eloss: DataBase + + def __init__(self, MT, xs, xs_offset, reference_frame, eloss): + type_ = REACTION_ELECTRON_EXCITATION + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.eloss = eloss + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + base = h5_group["energy_loss"] + eloss = DataTable(base["energy"][()], base["value"][()]) + + return cls(MT, xs, xs_offset, reference_frame, eloss) + + def __repr__(self): + text = super().__repr__() + text += f" - Energy loss: DataTable [ID: {self.eloss.ID}]\n" return text @@ -527,4 +620,4 @@ def set_energy_distribution(h5_group): else: print_error(f"Unsupported energy spectrum of type {spectrum_type}") - return energy_spectrum + return energy_spectrum \ No newline at end of file diff --git a/mcdc/object_/util.py b/mcdc/object_/util.py index dd63685b..ead6a10a 100644 --- a/mcdc/object_/util.py +++ b/mcdc/object_/util.py @@ -64,3 +64,465 @@ def is_sorted(a): # ====================================================================================== # Data # ====================================================================================== + +# Natural isotopic abundance data from +# https://www.nndc.bnl.gov/walletcards/search.html + +ISOTOPIC_ABUNDANCE = { + "H": { + "H1": (99.972 + 99.999) / 2, + "H2": (0.001 + 0.028) / 2, + }, + "He": { + "He3": 0.0002, + "He4": 99.9998, + }, + "Li": { + "Li6": (1.9 + 7.8) / 2, + "Li7": (92.2 + 98.1) / 2, + }, + "Be": { + "Be9": 100.0, + }, + "B": { + "B10": (18.9 + 20.4) / 2, + "B11": (79.6 + 81.1) / 2, + }, + "C": { + "C12": (98.84 + 99.04) / 2, + "C13": (0.96 + 1.16) / 2, + }, + "N": { + "N14": (99.578 + 99.663) / 2, + "N15": (0.337 + 0.422) / 2, + }, + "O": { + "O16": (99.738 + 99.776) / 2, + "O17": (0.0367 + 0.04) / 2, + "O18": (0.187 + 0.222) / 2, + }, + "F": { + "F19": 100.0, + }, + "Ne": { + "Ne20": 90.48, + "Ne21": 0.27, + "Ne22": 9.25, + }, + "Na": { + "Na23": 100.0, + }, + "Mg": { + "Mg24": (78.88 + 79.05) / 2, + "Mg25": (9.988 + 10.034) / 2, + "Mg26": (10.96 + 11.09) / 2, + }, + "Al": { + "Al27": 100.0, + }, + "Si": { + "Si28": (92.191 + 92.318) / 2, + "Si29": (4.645 + 4.699) / 2, + "Si30": (3.037 + 3.11) / 2, + }, + "P": { + "P31": 100.0, + }, + "S": { + "S32": (94.41 + 95.29) / 2, + "S33": (0.729 + 0.797) / 2, + "S34": (3.96 + 4.77) / 2, + "S36": (0.0129 + 0.0187) / 2, + }, + "Cl": { + "Cl35": (75.5 + 76.1) / 2, + "Cl37": (23.9 + 24.5) / 2, + }, + "Ar": { + "Ar36": 0.3336, + "Ar38": 0.0629, + "Ar40": 99.6035, + }, + "K": { + "K39": 93.2581, + "K40": 0.0117, + "K41": 6.7302, + }, + "Ca": { + "Ca40": 96.941, + "Ca42": 0.647, + "Ca43": 0.135, + "Ca44": 2.086, + "Ca46": 0.004, + "Ca48": 0.187, + }, + "Sc": { + "Sc45": 100.0, + }, + "Ti": { + "Ti46": 8.25, + "Ti47": 7.44, + "Ti48": 73.72, + "Ti49": 5.41, + "Ti50": 5.18, + }, + "V": { + "V50": 0.25, + "V51": 99.75, + }, + "Cr": { + "Cr50": 4.345, + "Cr52": 83.789, + "Cr53": 9.501, + "Cr54": 2.365, + }, + "Mn": { + "Mn55": 100.0, + }, + "Fe": { + "Fe54": 5.845, + "Fe56": 91.754, + "Fe57": 2.119, + "Fe58": 0.282, + }, + "Co": { + "Co59": 100.0, + }, + "Ni": { + "Ni58": 68.0769, + "Ni60": 26.2231, + "Ni61": 1.1399, + "Ni62": 3.6345, + "Ni64": 0.9256, + }, + "Cu": { + "Cu63": 69.15, + "Cu65": 30.85, + }, + "Zn": { + "Zn64": 49.17, + "Zn66": 27.73, + "Zn67": 4.04, + "Zn68": 18.45, + "Zn70": 0.61, + }, + "Ga": { + "Ga69": 60.108, + "Ga71": 39.892, + }, + "Ge": { + "Ge70": 20.52, + "Ge72": 27.45, + "Ge73": 7.76, + "Ge74": 36.52, + "Ge76": 7.75, + }, + "As": { + "As75": 100.0, + }, + "Se": { + "Se74": 0.86, + "Se76": 9.23, + "Se77": 7.6, + "Se78": 23.69, + "Se80": 49.8, + "Se82": 8.82, + }, + "Br": { + "Br79": (50.5 + 50.8) / 2, + "Br81": (49.2 + 49.5) / 2, + }, + "Kr": { + "Kr78": 0.355, + "Kr80": 2.286, + "Kr82": 11.593, + "Kr83": 11.5, + "Kr84": 56.987, + "Kr86": 17.279, + }, + "Rb": { + "Rb85": 72.17, + "Rb87": 27.83, + }, + "Sr": { + "Sr84": 0.56, + "Sr86": 9.86, + "Sr87": 7.0, + "Sr88": 82.58, + }, + "Y": { + "Y89": 100.0, + }, + "Zr": { + "Zr90": 51.45, + "Zr91": 11.22, + "Zr92": 17.15, + "Zr94": 17.38, + "Zr96": 2.8, + }, + "Nb": { + "Nb93": 100.0, + }, + "Mo": { + "Mo92": 14.649, + "Mo94": 9.187, + "Mo95": 15.873, + "Mo96": 16.673, + "Mo97": 9.582, + "Mo98": 24.292, + "Mo100": 9.744, + }, + "Ru": { + "Ru96": 5.54, + "Ru98": 1.87, + "Ru99": 12.76, + "Ru100": 12.6, + "Ru101": 17.06, + "Ru102": 31.55, + "Ru104": 18.62, + }, + "Rh": { + "Rh103": 100.0, + }, + "Pd": { + "Pd102": 1.02, + "Pd104": 11.14, + "Pd105": 22.33, + "Pd106": 27.33, + "Pd108": 26.46, + "Pd110": 11.72, + }, + "Ag": { + "Ag107": 51.839, + "Ag109": 48.161, + }, + "Cd": { + "Cd106": 1.245, + "Cd108": 0.888, + "Cd110": 12.47, + "Cd111": 12.795, + "Cd112": 24.109, + "Cd113": 12.227, + "Cd114": 28.754, + "Cd116": 7.512, + }, + "In": { + "In113": 4.281, + "In115": 95.719, + }, + "Sn": { + "Sn112": 0.97, + "Sn114": 0.66, + "Sn115": 0.34, + "Sn116": 14.54, + "Sn117": 7.68, + "Sn118": 24.22, + "Sn119": 8.59, + "Sn120": 32.58, + "Sn122": 4.63, + "Sn124": 5.79, + }, + "Sb": { + "Sb121": 57.21, + "Sb123": 42.79, + }, + "Te": { + "Te120": 0.09, + "Te122": 2.55, + "Te123": 0.89, + "Te124": 4.74, + "Te125": 7.07, + "Te126": 8.84, + "Te128": 31.74, + "Te130": 34.08, + }, + "I": { + "I127": 100.0, + }, + "Xe": { + "Xe124": 0.095, + "Xe126": 0.089, + "Xe128": 1.91, + "Xe129": 26.4, + "Xe130": 4.071, + "Xe131": 21.232, + "Xe132": 26.909, + "Xe134": 10.436, + "Xe136": 8.857, + }, + "Cs": { + "Cs133": 100.0, + }, + "Ba": { + "Ba130": 0.11, + "Ba132": 0.1, + "Ba134": 2.42, + "Ba135": 6.59, + "Ba136": 7.85, + "Ba137": 11.23, + "Ba138": 71.7, + }, + "La": { + "La138": 0.08881, + "La139": 99.91119, + }, + "Ce": { + "Ce136": 0.186, + "Ce138": 0.251, + "Ce140": 88.449, + "Ce142": 11.114, + }, + "Pr": { + "Pr141": 100.0, + }, + "Nd": { + "Nd142": 27.153, + "Nd143": 12.173, + "Nd144": 23.798, + "Nd145": 8.293, + "Nd146": 17.189, + "Nd148": 5.756, + "Nd150": 5.638, + }, + "Sm": { + "Sm144": 3.08, + "Sm147": 15.0, + "Sm149": 13.82, + "Sm150": 7.37, + "Sm152": 26.74, + "Sm154": 22.74, + }, + "Eu": { + "Eu151": 47.81, + "Eu153": 52.19, + }, + "Gd": { + "Gd152": 0.2, + "Gd154": 2.18, + "Gd155": 14.8, + "Gd156": 20.47, + "Gd157": 15.65, + "Gd158": 24.84, + "Gd160": 21.86, + }, + "Tb": { + "Tb159": 100.0, + }, + "Dy": { + "Dy156": 0.056, + "Dy158": 0.095, + "Dy160": 2.329, + "Dy161": 18.889, + "Dy162": 25.475, + "Dy163": 24.896, + "Dy164": 28.26, + }, + "Ho": { + "Ho165": 100.0, + }, + "Er": { + "Er162": 0.139, + "Er164": 1.601, + "Er166": 33.503, + "Er167": 22.869, + "Er168": 26.978, + "Er170": 14.91, + }, + "Tm": { + "Tm169": 100.0, + }, + "Yb": { + "Yb168": 0.123, + "Yb170": 2.982, + "Yb171": 14.086, + "Yb172": 21.686, + "Yb173": 16.103, + "Yb174": 32.025, + "Yb176": 12.995, + }, + "Lu": { + "Lu175": 97.401, + "Lu176": 2.599, + }, + "Hf": { + "Hf174": 0.16, + "Hf176": 5.26, + "Hf177": 18.6, + "Hf178": 27.28, + "Hf179": 13.62, + "Hf180": 35.08, + }, + "Ta": { + "Ta180": 0.01201, + "Ta181": 99.98799, + }, + "W": { + "W180": 0.12, + "W182": 26.5, + "W183": 14.31, + "W184": 30.64, + "W186": 28.43, + }, + "Re": { + "Re185": 37.4, + "Re187": 62.6, + }, + "Os": { + "Os184": 0.02, + "Os186": 1.59, + "Os187": 1.97, + "Os188": 13.24, + "Os189": 16.15, + "Os190": 26.26, + "Os192": 40.78, + }, + "Ir": { + "Ir191": 37.3, + "Ir193": 62.7, + }, + "Pt": { + "Pt190": 0.012, + "Pt192": 0.782, + "Pt194": 32.864, + "Pt195": 33.77, + "Pt196": 25.21, + "Pt198": 7.356, + }, + "Au": { + "Au197": 100.0, + }, + "Hg": { + "Hg196": 0.15, + "Hg198": 10.04, + "Hg199": 16.94, + "Hg200": 23.14, + "Hg201": 13.17, + "Hg202": 29.74, + "Hg204": 6.82, + }, + "Tl": { + "Tl203": (29.44 + 29.59) / 2, + "Tl205": (70.41 + 70.56) / 2, + }, + "Pb": { + "Pb204": 1.4, + "Pb206": 24.1, + "Pb207": 22.1, + "Pb208": 52.4, + }, + "Bi": { + "Bi209": 100.0, + }, + "Th": { + "Th230": 0.02, + "Th232": 99.98, + }, + "Pa": { + "Pa231": 100.0, + }, + "U": { + "U234": 0.0054, + "U235": 0.7204, + "U238": 99.2742, + }, +} diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py index 73da89d7..bcc4b284 100644 --- a/mcdc/transport/physics/electron/native.py +++ b/mcdc/transport/physics/electron/native.py @@ -1,6 +1,7 @@ import math import numpy as np +from mcdc.mcdc_get import reaction from numba import njit #### @@ -8,20 +9,16 @@ import mcdc.code_factory.adapt as adapt import mcdc.mcdc_get as mcdc_get import mcdc.numba_types as type_ -from mcdc.print_ import print_structure import mcdc.transport.particle as particle_module import mcdc.transport.particle_bank as particle_bank_module import mcdc.transport.rng as rng from mcdc.constant import ( - ANGLE_DISTRIBUTED, - ANGLE_ENERGY_CORRELATED, - ANGLE_ISOTROPIC, LIGHT_SPEED, PI, - PI_HALF, - PI_SQRT, + FINE_STRUCTURE_CONSTANT, ELECTRON_MASS, + ELECTRON_CUTOFF_ENERGY, REACTION_TOTAL, REACTION_ELECTRON_IONIZATION, REACTION_ELECTRON_ELASTIC_SCATTERING, @@ -30,14 +27,11 @@ ) from mcdc.transport.data import evaluate_data from mcdc.transport.distribution import ( - sample_correlated_distribution, sample_distribution, - sample_isotropic_cosine, - sample_isotropic_direction, sample_multi_table, ) from mcdc.transport.physics.util import evaluate_xs_energy_grid, scatter_direction -from mcdc.transport.util import find_bin, linear_interpolation +from mcdc.transport.util import linear_interpolation # ====================================================================================== @@ -65,14 +59,12 @@ def macro_xs(reaction_type, particle_container, mcdc, data): E = particle["E"] total = 0.0 - for i in range(material["N_element"]): element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) element = mcdc["elements"][element_ID] element_density = mcdc_get.native_material.element_densities(i, material, data) xs = total_micro_xs(reaction_type, E, element, data) - total += element_density * xs return total @@ -107,8 +99,7 @@ def reaction_micro_xs(E, reaction_base, element, data): offset = reaction_base["xs_offset_"] if idx < offset: return 0.0 - else: - idx -= offset + idx -= offset xs0 = mcdc_get.reaction.xs(idx, reaction_base, data) xs1 = mcdc_get.reaction.xs(idx + 1, reaction_base, data) @@ -124,14 +115,15 @@ def electron_production_xs(reaction_type, particle_container, mcdc, data): if reaction_type == REACTION_TOTAL: ionization_type = REACTION_ELECTRON_IONIZATION elastic_type = REACTION_ELECTRON_ELASTIC_SCATTERING - bremmstrahlung_type = REACTION_ELECTRON_BREMSSTRAHLUNG + bremsstrahlung_type = REACTION_ELECTRON_BREMSSTRAHLUNG excitation_type = REACTION_ELECTRON_EXCITATION ionization_xs = electron_production_xs(ionization_type, particle_container, mcdc, data) elastic_xs = electron_production_xs(elastic_type, particle_container, mcdc, data) - bremmstrahlung_xs = electron_production_xs(bremmstrahlung_type, particle_container, mcdc, data) + bremsstrahlung_xs = electron_production_xs(bremsstrahlung_type, particle_container, mcdc, data) excitation_xs = electron_production_xs(excitation_type, particle_container, mcdc, data) - return ionization_xs + elastic_xs + bremmstrahlung_xs + excitation_xs + + return ionization_xs + elastic_xs + bremsstrahlung_xs + excitation_xs elif reaction_type == REACTION_ELECTRON_IONIZATION: return 2 * macro_xs(reaction_type, particle_container, mcdc, data) @@ -152,3 +144,449 @@ def electron_production_xs(reaction_type, particle_container, mcdc, data): # ====================================================================================== # Collision # ====================================================================================== + + +@njit +def collision(particle_container, prog, data): + mcdc = adapt.mcdc_global(prog) + + particle = particle_container[0] + material = mcdc["native_materials"][particle["material_ID"]] + + # Particle properties + E = particle["E"] + + # Check for cutoff energy + if E <= ELECTRON_CUTOFF_ENERGY: + particle["alive"] = False + edep = E + return edep + + # ================================================================================== + # Sample colliding element + # ================================================================================== + + SigmaT = macro_xs(REACTION_TOTAL, particle_container, mcdc, data) + + xi = rng.lcg(particle_container) * SigmaT + total = 0.0 + + # Loop over elements in the material + for i in range(material["N_element"]): + element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) + element = mcdc["elements"][element_ID] + + element_density = mcdc_get.native_material.element_densities(i, material, data) + sigmaT = total_micro_xs(REACTION_TOTAL, E, element, data) + + total += element_density * sigmaT + + if total > xi: + break + + # ================================================================================== + # Sample and perform reaction + # ================================================================================== + + sigmaT_element = total_micro_xs(REACTION_TOTAL, E, element, data) + + sigma_ionization = total_micro_xs( + REACTION_ELECTRON_IONIZATION, E, element, data + ) + sigma_elastic = total_micro_xs( + REACTION_ELECTRON_ELASTIC_SCATTERING, E, element, data + ) + sigma_bremsstrahlung = total_micro_xs( + REACTION_ELECTRON_BREMSSTRAHLUNG, E, element, data + ) + sigma_excitation = total_micro_xs( + REACTION_ELECTRON_EXCITATION, E, element, data + ) + + xi = rng.lcg(particle_container) * sigmaT_element + total = 0.0 + + # Ionization + total += sigma_ionization + if xi < total: + total -= sigma_ionization + for i in range(element["N_ionization_reaction"]): + reaction_ID = int( + mcdc_get.element.ionization_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_ionization_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return ionization(reaction, particle_container, element, prog, data) + + # Elastic scattering + total += sigma_elastic + if xi < total: + total -= sigma_elastic + for i in range(element["N_elastic_scattering_reaction"]): + reaction_ID = int( + mcdc_get.element.elastic_scattering_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_elastic_scattering_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return elastic_scattering(reaction, particle_container, element, prog, data) + + # Bremsstrahlung + total += sigma_bremsstrahlung + if xi < total: + total -= sigma_bremsstrahlung + for i in range(element["N_bremsstrahlung_reaction"]): + reaction_ID = int( + mcdc_get.element.bremsstrahlung_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_bremsstrahlung_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return bremsstrahlung(reaction, particle_container, element, prog, data) + + + # Excitation + total += sigma_excitation + if xi < total: + total -= sigma_excitation + + for i in range(element["N_excitation_reaction"]): + reaction_ID = int( + mcdc_get.element.excitation_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_excitation_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return excitation(reaction, particle_container, element, prog, data) + + +# ====================================================================================== +# Elastic scattering +# ====================================================================================== + + +@njit +def compute_scattering_eta(E, Z): + pc = math.sqrt(E * (E + 2*ELECTRON_MASS)) + beta = pc / (E + ELECTRON_MASS) + tau = E / ELECTRON_MASS + + r = (FINE_STRUCTURE_CONSTANT * ELECTRON_MASS) / (0.885 * pc) + z_sq = float(Z) ** (2.0 / 3.0) + bracket = 1.13 + 3.76 * ((FINE_STRUCTURE_CONSTANT * float(Z)) / beta) ** 2 + rel = math.sqrt(tau / (tau + 1.0)) + + return 0.25 * (r * r) * z_sq * bracket * rel + + +@njit +def sample_small_angle_mu_coulomb(E, Z, rng_state, mu_cut): + eta = compute_scattering_eta(E, Z) + + x_cut = 1.0 - mu_cut + u = rng.lcg(rng_state) + + denom = (1.0 / eta) - (1.0 / (eta + x_cut)) + inv = (1.0 / eta) - u * denom + x = (1.0 / inv) - eta + + return 1.0 - x + + +@njit +def elastic_scattering(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + # ------------------------------------------------------------------------- + # Total elastic xs + # ------------------------------------------------------------------------- + reaction_base_ID = int(reaction["parent_ID"]) + reaction_base = mcdc["reactions"][reaction_base_ID] + xs_total = reaction_micro_xs(E, reaction_base, element, data) + + # If large-angle, xs from data table (EEDL MT=525) + xs_large = 0.0 + xs_large_ID = int(reaction["xs_large_ID"]) + + if xs_large_ID >= 0: + xs_large = elastic_large_xs(E, reaction, mcdc, data) + + # Important to check because of numerical issues + if xs_large < 0.0: + xs_large = 0.0 + if xs_large > xs_total: + xs_large = xs_total + + prob_large = xs_large / xs_total + mu_cut = float(reaction["mu_cut"]) + + xi = rng.lcg(particle_container) + + if xi < prob_large: + # --------------------------------------------------------------------- + # Large-angle elastic scattering + # --------------------------------------------------------------------- + + # Sample mu from EEDL tabulated distribution + mu_table_ID = int(reaction["mu_table_ID"]) + multi_table = mcdc["multi_table_distributions"][mu_table_ID] + mu0 = sample_multi_table(E, particle_container, multi_table, data) + + else: + # --------------------------------------------------------------------- + # Small-angle elastic scattering (Coulomb tail sampling) + # --------------------------------------------------------------------- + + Z = int(element["atomic_number"]) + mu0 = sample_small_angle_mu_coulomb(E, Z, particle_container, mu_cut) + + # Update direction + azi = 2.0 * PI * rng.lcg(particle_container) + # Current direction + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) + + particle["ux"] = ux_new + particle["uy"] = uy_new + particle["uz"] = uz_new + + return 0.0 + + +@njit +def elastic_large_xs(E, reaction, mcdc, data): + xs_id = int(reaction["xs_large_ID"]) + data_base = mcdc["data"][xs_id] + return evaluate_data(E, data_base, mcdc, data) + + +# ============================================================================== +# Excitation +# ============================================================================== + + +@njit +def excitation(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + dE = evaluate_eloss(E, reaction, mcdc, data) + + # Calculate outgoing energy + E_out = E - dE + + # Check for cutoff + if E_out <= ELECTRON_CUTOFF_ENERGY: + # deposit remaining and kill + edep = E + particle["E"] = 0.0 + particle["alive"] = False + return edep + + # If above cutoff, just deposit dE + particle["E"] = E_out + edep = dE + return edep + + +@njit +def evaluate_eloss(E, reaction, mcdc, data): + data_base = mcdc["data"][int(reaction["eloss_ID"])] + return evaluate_data(E, data_base, mcdc, data) + + +# ============================================================================== +# Bremsstrahlung (photon not tracked - energy deposit = 0 ) +# ============================================================================== + + +@njit +def bremsstrahlung(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + dE = evaluate_eloss(E, reaction, mcdc, data) + + E_out = E - dE + # Check for cutoff + if E_out <= ELECTRON_CUTOFF_ENERGY: + # deposit remaining and kill + edep = E_out + particle["E"] = 0.0 + particle["alive"] = False + return edep + + # If above cutoff, allow photon to escape, just update electron energy + particle["E"] = E_out + edep = 0.0 + return edep + + +# ============================================================================== +# Ionization +# ============================================================================== + + +@njit +def ionization(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + # Sample subshell + N = int(reaction["N_subshell"]) + xs_vals = np.empty(N, dtype=np.float64) + + total = 0.0 + for i in range(N): + xs_sub_id = int(reaction["subshell_xs_IDs"][i]) + xs_sub_table = mcdc["data"][xs_sub_id] + xs_sub_i = evaluate_data(E, xs_sub_table, mcdc, data) + xs_vals[i] = xs_sub_i + total += xs_sub_i + + xi = rng.lcg(particle_container) * total + total_acc = 0.0 + chosen = 0 + for i in range(N): + total_acc += xs_vals[i] + if total_acc >= xi: + chosen = i + break + + # Binding energy (from element, shared across reactions) + B = mcdc_get.element.ionization_subshell_binding_energy(chosen, element, data) + # Deposit all energy if below binding energy + if E <= B: + edep = E + particle["alive"] = False + particle["E"] = 0.0 + return edep + + # Sample secondary energy T_delta from distribution + dist_ID = int(reaction["subshell_product_IDs"][chosen]) + dist_base = mcdc["distributions"][dist_ID] + T_delta = sample_distribution(E, dist_base, particle_container, mcdc, data, scale=True) + + # Primary outgoing energy + E_out = E - B - T_delta + particle["E"] = E_out + + edep = B + + primary_alive_after = True + if E_out <= ELECTRON_CUTOFF_ENERGY: + edep += E_out + particle["E"] = 0.0 + particle["alive"] = False + primary_alive_after = False + + if T_delta <= ELECTRON_CUTOFF_ENERGY: + edep += T_delta + return edep + + # Sample delta direction + ux_delta, uy_delta, uz_delta = sample_delta_direction(T_delta, E, particle_container) + + # Momentum conservation if primary alive + if primary_alive_after: + p_before = math.sqrt(E * (E + 2*ELECTRON_MASS)) + p_delta = math.sqrt(T_delta * (T_delta + 2*ELECTRON_MASS)) + + ux_before = particle["ux"] + uy_before = particle["uy"] + uz_before = particle["uz"] + + # Momentum vectors after collision + px_after = p_before * ux_before - p_delta * ux_delta + py_after = p_before * uy_before - p_delta * uy_delta + pz_after = p_before * uz_before - p_delta * uz_delta + + # Normalize and set primary's new direction + norm_sq = px_after * px_after + py_after * py_after + pz_after * pz_after + if norm_sq > 0.0: + norm = math.sqrt(norm_sq) + particle["ux"] = px_after / norm + particle["uy"] = py_after / norm + particle["uz"] = pz_after / norm + + # Add secondary particle to bank + particle_container_new = np.zeros(1, type_.particle_data) + particle_new = particle_container_new[0] + particle_module.copy_as_child(particle_container_new, particle_container) + + particle_new["E"] = T_delta + particle_new["ux"] = ux_delta + particle_new["uy"] = uy_delta + particle_new["uz"] = uz_delta + particle_new["w"] = particle["w"] + + particle_bank_module.add_active(particle_container_new, prog) + return edep + + +@njit +def compute_mu_delta(T_delta, T_prim): + me = ELECTRON_MASS + pd = math.sqrt(T_delta * (T_delta + 2.0 * me)) + pp = math.sqrt(T_prim * (T_prim + 2.0 * me)) + mu = (T_delta * (T_prim + 2.0 * me)) / (pd * pp) + # Check in case of numerical issues + if mu < -1.0: + mu = -1.0 + if mu > 1.0: + mu = 1.0 + + return mu + + +@njit +def sample_delta_direction(T_delta, T_prim, particle_container): + mu = compute_mu_delta(T_delta, T_prim) + phi = 2.0 * PI * rng.lcg(particle_container) + s = math.sqrt(1.0 - mu * mu) + ux = s * math.cos(phi) + uy = s * math.sin(phi) + uz = mu + return ux, uy, uz + + +@njit +def subshell_xs(E, xs_id, mcdc, data): + data_base = mcdc["data"][int(xs_id)] + return evaluate_data(E, data_base, mcdc, data) From 426c5164ece27beaa0a8e52767556e93664ef733 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:41:56 +0700 Subject: [PATCH 09/22] add electron to physics __init__ --- mcdc/transport/physics/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mcdc/transport/physics/__init__.py b/mcdc/transport/physics/__init__.py index 73dc04af..66133009 100644 --- a/mcdc/transport/physics/__init__.py +++ b/mcdc/transport/physics/__init__.py @@ -5,4 +5,5 @@ collision_distance, collision, ) +import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron From fd995dc0eea21a0582eb61d5bc6201d1ac0abc7b Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:42:25 +0700 Subject: [PATCH 10/22] add electron to physics interface --- mcdc/transport/physics/interface.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py index fd30f0f8..e9df90c9 100644 --- a/mcdc/transport/physics/interface.py +++ b/mcdc/transport/physics/interface.py @@ -5,6 +5,7 @@ ### import mcdc.transport.rng as rng +import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron from mcdc.constant import * @@ -20,6 +21,8 @@ def particle_speed(particle_container, mcdc, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.particle_speed(particle_container, mcdc, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + return electron.particle_speed(particle_container, mcdc, data) return -1.0 @@ -33,17 +36,22 @@ def macro_xs(reaction_type, particle_container, mcdc, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.macro_xs(reaction_type, particle_container, mcdc, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + return electron.macro_xs(reaction_type, particle_container, mcdc, data) return -1.0 @njit def neutron_production_xs(reaction_type, particle_container, mcdc, data): + """ + This is currently only used for neutron k-eigenvalue simulation + """ particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.neutron_production_xs( reaction_type, particle_container, mcdc, data ) - return -1.0 + return 0.0 # ====================================================================================== @@ -71,3 +79,5 @@ def collision(particle_container, prog, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: neutron.collision(particle_container, prog, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + electron.collision(particle_container, prog, data) From eb4c4fb3f4e6c8b5a123403ceb860b1c2e0c0795 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:42:57 +0700 Subject: [PATCH 11/22] minor edits for neutron physics --- mcdc/transport/physics/neutron/native.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/mcdc/transport/physics/neutron/native.py b/mcdc/transport/physics/neutron/native.py index 99f6f3c4..433ff9da 100644 --- a/mcdc/transport/physics/neutron/native.py +++ b/mcdc/transport/physics/neutron/native.py @@ -237,6 +237,7 @@ def collision(particle_container, prog, data): particle["w"] *= (SigmaT - SigmaC) / SigmaT SigmaT -= SigmaC + # Loop over nuclides in the material xi = rng.lcg(particle_container) * SigmaT total = 0.0 for i in range(material["N_nuclide"]): @@ -251,8 +252,7 @@ def collision(particle_container, prog, data): particle["w"] *= (sigmaT - sigmaC) / sigmaT sigmaT -= sigmaC - SigmaT_nuclide = nuclide_density * sigmaT - total += SigmaT_nuclide + total += nuclide_density * sigmaT if total > xi: break @@ -260,6 +260,9 @@ def collision(particle_container, prog, data): # ================================================================================== # Sample and perform reaction # ================================================================================== + # TODO: The switches for the different reactions have a common code pattern. + # Modularizing the pattern would improve maintainability. + # Note that this also applies to other particle physics. sigma_elastic = total_micro_xs( REACTION_NEUTRON_ELASTIC_SCATTERING, E, nuclide, data @@ -299,7 +302,6 @@ def collision(particle_container, prog, data): total += sigma_inelastic if xi < total: total -= sigma_inelastic - for i in range(nuclide["N_inelastic_scattering_reaction"]): reaction_ID = int( mcdc_get.nuclide.inelastic_scattering_reaction_IDs(i, nuclide, data) @@ -307,8 +309,7 @@ def collision(particle_container, prog, data): reaction = mcdc["neutron_inelastic_scattering_reactions"][reaction_ID] reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] - xs = reaction_micro_xs(E, reaction_base, nuclide, data) - total += xs + total += reaction_micro_xs(E, reaction_base, nuclide, data) if xi < total: inelastic_scattering(reaction, particle_container, nuclide, prog, data) return From 0ad930921bf800213445961b259eaf1527a17633 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:43:39 +0700 Subject: [PATCH 12/22] correct electron physics __init__ --- mcdc/transport/physics/electron/__init__.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/mcdc/transport/physics/electron/__init__.py b/mcdc/transport/physics/electron/__init__.py index 4715cd68..7a0a5f2d 100644 --- a/mcdc/transport/physics/electron/__init__.py +++ b/mcdc/transport/physics/electron/__init__.py @@ -1,3 +1,6 @@ -import mcdc.physics.electron.interface as interface - -import mcdc.physics.electron.native as native \ No newline at end of file +from .interface import ( + particle_speed, + macro_xs, + collision, +) +import mcdc.transport.physics.electron.native as native From 846423399bbb01b2e80f05d1542ad62ecf21fedc Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:44:02 +0700 Subject: [PATCH 13/22] correct electron physics interface. Also delete electron_production_xs --- mcdc/transport/physics/electron/interface.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/mcdc/transport/physics/electron/interface.py b/mcdc/transport/physics/electron/interface.py index f489dc1a..8fc991ff 100644 --- a/mcdc/transport/physics/electron/interface.py +++ b/mcdc/transport/physics/electron/interface.py @@ -11,7 +11,7 @@ @njit -def particle_speed(particle_container, material, data): +def particle_speed(particle_container, mcdc, data): return native.particle_speed(particle_container) @@ -21,13 +21,8 @@ def particle_speed(particle_container, material, data): @njit -def macro_xs(reaction_type, material, particle_container, mcdc, data): - return native.macro_xs(reaction_type, material, particle_container, mcdc, data) - - -@njit -def electron_production_xs(reaction_type, particle_container, mcdc, data): - return native.electron_production_xs(reaction_type, particle_container, mcdc, data) +def macro_xs(reaction_type, particle_container, mcdc, data): + return native.macro_xs(reaction_type, particle_container, mcdc, data) # ====================================================================================== @@ -37,4 +32,4 @@ def electron_production_xs(reaction_type, particle_container, mcdc, data): @njit def collision(particle_container, prog, data): - return native.collision(particle_container, prog, data) \ No newline at end of file + return native.collision(particle_container, prog, data) From 90110fc6f15b36f215e104833b5cf0ea2751018c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:45:32 +0700 Subject: [PATCH 14/22] minor edits on electron native physics. --- mcdc/transport/physics/electron/native.py | 102 ++++++++-------------- 1 file changed, 34 insertions(+), 68 deletions(-) diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py index bcc4b284..ee740292 100644 --- a/mcdc/transport/physics/electron/native.py +++ b/mcdc/transport/physics/electron/native.py @@ -106,41 +106,6 @@ def reaction_micro_xs(E, reaction_base, element, data): return linear_interpolation(E, E0, E1, xs0, xs1) -@njit -def electron_production_xs(reaction_type, particle_container, mcdc, data): - particle = particle_container[0] - material_base = mcdc["materials"][particle["material_ID"]] - material = mcdc["native_materials"][material_base["child_ID"]] - - if reaction_type == REACTION_TOTAL: - ionization_type = REACTION_ELECTRON_IONIZATION - elastic_type = REACTION_ELECTRON_ELASTIC_SCATTERING - bremsstrahlung_type = REACTION_ELECTRON_BREMSSTRAHLUNG - excitation_type = REACTION_ELECTRON_EXCITATION - - ionization_xs = electron_production_xs(ionization_type, particle_container, mcdc, data) - elastic_xs = electron_production_xs(elastic_type, particle_container, mcdc, data) - bremsstrahlung_xs = electron_production_xs(bremsstrahlung_type, particle_container, mcdc, data) - excitation_xs = electron_production_xs(excitation_type, particle_container, mcdc, data) - - return ionization_xs + elastic_xs + bremsstrahlung_xs + excitation_xs - - elif reaction_type == REACTION_ELECTRON_IONIZATION: - return 2 * macro_xs(reaction_type, particle_container, mcdc, data) - - elif reaction_type == REACTION_ELECTRON_ELASTIC_SCATTERING: - return macro_xs(reaction_type, particle_container, mcdc, data) - - elif reaction_type == REACTION_ELECTRON_BREMSSTRAHLUNG: - return macro_xs(reaction_type, particle_container, mcdc, data) - - elif reaction_type == REACTION_ELECTRON_EXCITATION: - return macro_xs(reaction_type, particle_container, mcdc, data) - - else: - return -1.0 - - # ====================================================================================== # Collision # ====================================================================================== @@ -168,10 +133,9 @@ def collision(particle_container, prog, data): SigmaT = macro_xs(REACTION_TOTAL, particle_container, mcdc, data) + # Loop over elements in the material xi = rng.lcg(particle_container) * SigmaT total = 0.0 - - # Loop over elements in the material for i in range(material["N_element"]): element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) element = mcdc["elements"][element_ID] @@ -187,23 +151,20 @@ def collision(particle_container, prog, data): # ================================================================================== # Sample and perform reaction # ================================================================================== + # TODO: The switches for the different reactions have a common code pattern. + # Modularizing the pattern would improve maintainability. + # Note that this also applies to other particle physics. - sigmaT_element = total_micro_xs(REACTION_TOTAL, E, element, data) - - sigma_ionization = total_micro_xs( - REACTION_ELECTRON_IONIZATION, E, element, data - ) + sigma_ionization = total_micro_xs(REACTION_ELECTRON_IONIZATION, E, element, data) sigma_elastic = total_micro_xs( REACTION_ELECTRON_ELASTIC_SCATTERING, E, element, data ) sigma_bremsstrahlung = total_micro_xs( REACTION_ELECTRON_BREMSSTRAHLUNG, E, element, data - ) - sigma_excitation = total_micro_xs( - REACTION_ELECTRON_EXCITATION, E, element, data - ) + ) + sigma_excitation = total_micro_xs(REACTION_ELECTRON_EXCITATION, E, element, data) - xi = rng.lcg(particle_container) * sigmaT_element + xi = rng.lcg(particle_container) * sigmaT total = 0.0 # Ionization @@ -218,10 +179,10 @@ def collision(particle_container, prog, data): reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] total += reaction_micro_xs(E, reaction_base, element, data) - + if xi < total: return ionization(reaction, particle_container, element, prog, data) - + # Elastic scattering total += sigma_elastic if xi < total: @@ -234,10 +195,12 @@ def collision(particle_container, prog, data): reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] total += reaction_micro_xs(E, reaction_base, element, data) - + if xi < total: - return elastic_scattering(reaction, particle_container, element, prog, data) - + return elastic_scattering( + reaction, particle_container, element, prog, data + ) + # Bremsstrahlung total += sigma_bremsstrahlung if xi < total: @@ -250,11 +213,10 @@ def collision(particle_container, prog, data): reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] total += reaction_micro_xs(E, reaction_base, element, data) - + if xi < total: return bremsstrahlung(reaction, particle_container, element, prog, data) - # Excitation total += sigma_excitation if xi < total: @@ -268,7 +230,7 @@ def collision(particle_container, prog, data): reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] total += reaction_micro_xs(E, reaction_base, element, data) - + if xi < total: return excitation(reaction, particle_container, element, prog, data) @@ -280,8 +242,8 @@ def collision(particle_container, prog, data): @njit def compute_scattering_eta(E, Z): - pc = math.sqrt(E * (E + 2*ELECTRON_MASS)) - beta = pc / (E + ELECTRON_MASS) + pc = math.sqrt(E * (E + 2 * ELECTRON_MASS)) + beta = pc / (E + ELECTRON_MASS) tau = E / ELECTRON_MASS r = (FINE_STRUCTURE_CONSTANT * ELECTRON_MASS) / (0.885 * pc) @@ -317,7 +279,7 @@ def elastic_scattering(reaction, particle_container, element, prog, data): E = particle["E"] # ------------------------------------------------------------------------- - # Total elastic xs + # Total elastic xs # ------------------------------------------------------------------------- reaction_base_ID = int(reaction["parent_ID"]) reaction_base = mcdc["reactions"][reaction_base_ID] @@ -467,7 +429,7 @@ def ionization(reaction, particle_container, element, prog, data): # Current energy E = particle["E"] - + # Sample subshell N = int(reaction["N_subshell"]) xs_vals = np.empty(N, dtype=np.float64) @@ -497,11 +459,13 @@ def ionization(reaction, particle_container, element, prog, data): particle["alive"] = False particle["E"] = 0.0 return edep - + # Sample secondary energy T_delta from distribution dist_ID = int(reaction["subshell_product_IDs"][chosen]) dist_base = mcdc["distributions"][dist_ID] - T_delta = sample_distribution(E, dist_base, particle_container, mcdc, data, scale=True) + T_delta = sample_distribution( + E, dist_base, particle_container, mcdc, data, scale=True + ) # Primary outgoing energy E_out = E - B - T_delta @@ -519,14 +483,16 @@ def ionization(reaction, particle_container, element, prog, data): if T_delta <= ELECTRON_CUTOFF_ENERGY: edep += T_delta return edep - + # Sample delta direction - ux_delta, uy_delta, uz_delta = sample_delta_direction(T_delta, E, particle_container) + ux_delta, uy_delta, uz_delta = sample_delta_direction( + T_delta, E, particle_container + ) # Momentum conservation if primary alive if primary_alive_after: - p_before = math.sqrt(E * (E + 2*ELECTRON_MASS)) - p_delta = math.sqrt(T_delta * (T_delta + 2*ELECTRON_MASS)) + p_before = math.sqrt(E * (E + 2 * ELECTRON_MASS)) + p_delta = math.sqrt(T_delta * (T_delta + 2 * ELECTRON_MASS)) ux_before = particle["ux"] uy_before = particle["uy"] @@ -544,7 +510,7 @@ def ionization(reaction, particle_container, element, prog, data): particle["ux"] = px_after / norm particle["uy"] = py_after / norm particle["uz"] = pz_after / norm - + # Add secondary particle to bank particle_container_new = np.zeros(1, type_.particle_data) particle_new = particle_container_new[0] @@ -569,8 +535,8 @@ def compute_mu_delta(T_delta, T_prim): # Check in case of numerical issues if mu < -1.0: mu = -1.0 - if mu > 1.0: - mu = 1.0 + if mu > 1.0: + mu = 1.0 return mu From 41385263d15eac61e05b471f600d25439b943f47 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 18 Dec 2025 20:58:38 +0700 Subject: [PATCH 15/22] minor adjustments to excitation and bremsstrahlung --- mcdc/transport/physics/electron/native.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py index ee740292..b118872b 100644 --- a/mcdc/transport/physics/electron/native.py +++ b/mcdc/transport/physics/electron/native.py @@ -344,12 +344,12 @@ def elastic_large_xs(E, reaction, mcdc, data): # ============================================================================== -# Excitation +# Excitation (photon not tracked - all energy deposited) # ============================================================================== @njit -def excitation(reaction, particle_container, element, prog, data): +def excitation(reaction, particle_container, prog, data): mcdc = adapt.mcdc_global(prog) # Particle attributes @@ -389,7 +389,7 @@ def evaluate_eloss(E, reaction, mcdc, data): @njit -def bremsstrahlung(reaction, particle_container, element, prog, data): +def bremsstrahlung(reaction, particle_container, prog, data): mcdc = adapt.mcdc_global(prog) # Particle attributes From 14af77a4a4a11271d68218cf0afa5c8f39d32a4d Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Fri, 19 Dec 2025 15:52:31 -0800 Subject: [PATCH 16/22] add an example (needs to update) and refs --- .../lockwood/experimental_conf_Al.py | 133 ++++++++++++++++++ test/regression/electron_tests/references.md | 7 + 2 files changed, 140 insertions(+) create mode 100644 test/regression/electron_tests/lockwood/experimental_conf_Al.py create mode 100644 test/regression/electron_tests/references.md diff --git a/test/regression/electron_tests/lockwood/experimental_conf_Al.py b/test/regression/electron_tests/lockwood/experimental_conf_Al.py new file mode 100644 index 00000000..00663d6d --- /dev/null +++ b/test/regression/electron_tests/lockwood/experimental_conf_Al.py @@ -0,0 +1,133 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +# Set the XS library directory +os.environ["MCDC_XSLIB"] = os.path.join(os.path.dirname(os.path.dirname(os.getcwd())), "dummy_data") + +# ============================================================================= +# Set problem parameters +# ============================================================================= +MATERIAL_SYMBOL = "Al" +CSDA_RANGE = 0.569 # g/cm2 +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Other parameters (don't change) +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +TINY = 1e-12 +VOID_SCALE = 1.0e-30 +VOID_DENSITY = MAT_DENSITY_ATOMS_PER_BARN_CM * VOID_SCALE + +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) + +# ============================================================================= +# Exactly the same computational configuration as MCNP for Lockwood exp (from the MCNP VV slides) +# https://mcnp.lanl.gov/pdf_files/TechReport_2024_LANL_LA-UR-23-32743Rev.1_Kulesza.pdf +# ============================================================================= +# pz 0.0 +# pz 0.0233 - front foil thickness = 0.0233 cm +# pz 0.1233 - +0.1 cm void gap +# pz 0.12517037037 - calorimeter foil thickness = 0.00187037037 cm (5.05e-3/2.7) +# pz 0.22517037037 - +0.1 cm void gap +# pz 5.0 - thick Al plate end +# +# cz 100.0 - radius 100 cm boundary +z0 = 0.0 +z1 = 0.0233 +z2 = 0.1233 +z3 = 0.12517037037037037 +z4 = 0.22517037037037037 +z5 = 5.0 + +R_CYL = 100.0 + +# ============================================================================= +# Set materials +# ============================================================================= +mat_al = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +mat_void = mcdc.Material( + element_composition={MATERIAL_SYMBOL: VOID_DENSITY} +) + +# ============================================================================= +# Surfaces +# ============================================================================= +# Z-direction surfaces for layers + +s0 = mcdc.surface("plane-z", z=z0, bc="vacuum") +s1 = mcdc.surface("plane-z", z=z1) +s2 = mcdc.surface("plane-z", z=z2) +s3 = mcdc.surface("plane-z", z=z3) +s4 = mcdc.surface("plane-z", z=z4) +s5 = mcdc.surface("plane-z", z=z5, bc="vacuum") + +# ============================================================================= +# Cells +# ============================================================================= + +# Front foil (Al): z0 < z < z1 +mcdc.cell(+s0 & -s1, mat_al) + +# Void gap (almost vacuum): z1 < z < z2 +mcdc.cell(+s1 & -s2, mat_void) + +# Calorimeter foil (Al): z2 < z < z3 +mcdc.cell(+s2 & -s3, mat_al) + +# Void gap: z3 < z < z4 +mcdc.cell(+s3 & -s4, mat_void) + +# Thick Al plate: z4 < z < z5 +mcdc.cell(+s4 & -s5, mat_al) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0, 0 degree incidence +theta = 0.0 +mcdc.source( + z=[z0 + 1e-8, z0 + 1e-8], + particle_type="electron", + energy=np.array([[1e6 - 1.0, 1e6 + 1.0], [0.5, 0.5]]), + direction=[math.sin(theta), TINY, math.cos(theta)], +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(z4, z4+L, 101) + +mcdc.tally.mesh_tally( + scores=["edep", "flux"], + z=z_bins +) + +mcdc.tally.surface_tally(s0, scores=["net-current"]) +mcdc.tally.surface_tally(s5, scores=["net-current"]) + + +# ============================================================================= +# Settings and run +# ============================================================================= +settings = mcdc.Settings( + N_particle=100, + active_bank_buffer=10000 +) + +settings.save_input_deck = True +settings.output_name = f"lockwood_exp_setup_{datetime.now():%Y%m%d_%H%M%S}" +settings.debug_energy = False +settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/test/regression/electron_tests/references.md b/test/regression/electron_tests/references.md new file mode 100644 index 00000000..22e07924 --- /dev/null +++ b/test/regression/electron_tests/references.md @@ -0,0 +1,7 @@ +# References used during this work. + +1 - [Calorimetric Measurement of Electron Energy Deposition In Extended Media Theory vs Experiment - Lockwood](https://www.osti.gov/servlets/purl/5516784) +1 - [Validation of MCNP6 for Electron Energy Deposition in Extended Media](https://mcnp.lanl.gov/pdf_files/TechReport_2015_LANL_LA-UR-15-25134_DixonHughes.pdf) +2 - [An Electron/Photon/Relaxation Data Library for MCNP6](https://nucleardata.lanl.gov/files/la-ur-13-27377.pdf) +3 - [Electron/Photon/Relaxation Data Library—eprdata14](https://nucleardata.lanl.gov/ace/eprdata14/) +4 - [MCNP6.3 Electron Energy Deposition Validation with the Lockwood Experiments](https://mcnp.lanl.gov/pdf_files/TechReport_2024_LANL_LA-UR-23-32743Rev.1_Kulesza.pdf) \ No newline at end of file From f65e0d7821c162bb87b3757447a7e310f0bdff8d Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Fri, 19 Dec 2025 15:53:50 -0800 Subject: [PATCH 17/22] change elastic scattering helpers location --- mcdc/transport/physics/electron/native.py | 56 +++++++++++------------ 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py index b118872b..bbef71b5 100644 --- a/mcdc/transport/physics/electron/native.py +++ b/mcdc/transport/physics/electron/native.py @@ -240,34 +240,6 @@ def collision(particle_container, prog, data): # ====================================================================================== -@njit -def compute_scattering_eta(E, Z): - pc = math.sqrt(E * (E + 2 * ELECTRON_MASS)) - beta = pc / (E + ELECTRON_MASS) - tau = E / ELECTRON_MASS - - r = (FINE_STRUCTURE_CONSTANT * ELECTRON_MASS) / (0.885 * pc) - z_sq = float(Z) ** (2.0 / 3.0) - bracket = 1.13 + 3.76 * ((FINE_STRUCTURE_CONSTANT * float(Z)) / beta) ** 2 - rel = math.sqrt(tau / (tau + 1.0)) - - return 0.25 * (r * r) * z_sq * bracket * rel - - -@njit -def sample_small_angle_mu_coulomb(E, Z, rng_state, mu_cut): - eta = compute_scattering_eta(E, Z) - - x_cut = 1.0 - mu_cut - u = rng.lcg(rng_state) - - denom = (1.0 / eta) - (1.0 / (eta + x_cut)) - inv = (1.0 / eta) - u * denom - x = (1.0 / inv) - eta - - return 1.0 - x - - @njit def elastic_scattering(reaction, particle_container, element, prog, data): mcdc = adapt.mcdc_global(prog) @@ -336,6 +308,34 @@ def elastic_scattering(reaction, particle_container, element, prog, data): return 0.0 +@njit +def compute_scattering_eta(E, Z): + pc = math.sqrt(E * (E + 2 * ELECTRON_MASS)) + beta = pc / (E + ELECTRON_MASS) + tau = E / ELECTRON_MASS + + r = (FINE_STRUCTURE_CONSTANT * ELECTRON_MASS) / (0.885 * pc) + z_sq = float(Z) ** (2.0 / 3.0) + bracket = 1.13 + 3.76 * ((FINE_STRUCTURE_CONSTANT * float(Z)) / beta) ** 2 + rel = math.sqrt(tau / (tau + 1.0)) + + return 0.25 * (r * r) * z_sq * bracket * rel + + +@njit +def sample_small_angle_mu_coulomb(E, Z, rng_state, mu_cut): + eta = compute_scattering_eta(E, Z) + + x_cut = 1.0 - mu_cut + u = rng.lcg(rng_state) + + denom = (1.0 / eta) - (1.0 / (eta + x_cut)) + inv = (1.0 / eta) - u * denom + x = (1.0 / inv) - eta + + return 1.0 - x + + @njit def elastic_large_xs(E, reaction, mcdc, data): xs_id = int(reaction["xs_large_ID"]) From e362d69c84ff27b7c4a218fd92c42f18e2f49cae Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Fri, 19 Dec 2025 15:57:22 -0800 Subject: [PATCH 18/22] add slab example for Al --- .../electron_tests/lockwood/slab_Al.py | 101 ++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 test/regression/electron_tests/lockwood/slab_Al.py diff --git a/test/regression/electron_tests/lockwood/slab_Al.py b/test/regression/electron_tests/lockwood/slab_Al.py new file mode 100644 index 00000000..607af8dc --- /dev/null +++ b/test/regression/electron_tests/lockwood/slab_Al.py @@ -0,0 +1,101 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +# Set the XS library directory +os.environ["MCDC_XSLIB"] = os.path.join(os.path.dirname(os.path.dirname(os.getcwd())), "dummy_data") + +# ============================================================================= +# Set problem parameters +# ============================================================================= +MATERIAL_SYMBOL = "Al" +CSDA_RANGE = 0.569 # g/cm2 +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Other parameters (don't change) +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +TINY = 1e-12 +VOID_SCALE = 1.0e-30 +VOID_DENSITY = MAT_DENSITY_ATOMS_PER_BARN_CM * VOID_SCALE +M = 20 + +L = CSDA_RANGE / RHO_G_CM3 * M # cm +N_LAYERS = int(L / dz) + +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.surface("plane-z", z=0, bc="vacuum") +s2 = mcdc.surface("plane-z", z=L, bc="vacuum") + +mcdc.cell(+s1 & -s2, mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 +theta = 0 +z0 = 0.0 + +mcdc.source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[1e6 - 1, 1e6 + 1], [0.5, 0.5]]), + direction=[math.sin(theta), 0.0+TINY, math.cos(theta)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L/M, N_LAYERS + 1) + +mcdc.tally.mesh_tally( + scores=["edep", "flux"], + z=z_bins +) + +mcdc.tally.surface_tally(s1, scores=["net-current"]) +mcdc.tally.surface_tally(s2, scores=["net-current"]) + + +# ============================================================================= +# Settings and run +# ============================================================================= +settings = mcdc.Settings( + N_particle=100, + active_bank_buffer=10000 +) + +settings.save_input_deck = True +settings.output_name = f"lockwood_output_{datetime.now():%Y%m%d_%H%M%S}" +settings.debug_energy = False +settings.use_progress_bar = True + + +mcdc.run() \ No newline at end of file From a8eefaaaef3b4d1604c2fa74340d9b5602dc71ba Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Sat, 20 Dec 2025 16:08:47 -0800 Subject: [PATCH 19/22] examples that Ilham uploaded --- examples/kobayashi/input.py | 75 ++++++++++++++++++++++++++++ examples/kobayashi/process-output.py | 48 ++++++++++++++++++ 2 files changed, 123 insertions(+) create mode 100644 examples/kobayashi/input.py create mode 100644 examples/kobayashi/process-output.py diff --git a/examples/kobayashi/input.py b/examples/kobayashi/input.py new file mode 100644 index 00000000..8f959885 --- /dev/null +++ b/examples/kobayashi/input.py @@ -0,0 +1,75 @@ +import numpy as np +import mcdc + + +# ====================================================================================== +# Set model +# ====================================================================================== +# Based on Kobayashi dog-leg benchmark problem +# (PNE 2001, https://doi.org/10.1016/S0149-1970(01)00007-5) + +# Set materials +m = mcdc.MaterialMG(capture=np.array([0.05]), scatter=np.array([[0.05]])) +m_void = mcdc.MaterialMG(capture=np.array([5e-5]), scatter=np.array([[5e-5]])) + +# Set surfaces +sx1 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="reflective") +sx2 = mcdc.Surface.PlaneX(x=10.0) +sx3 = mcdc.Surface.PlaneX(x=30.0) +sx4 = mcdc.Surface.PlaneX(x=40.0) +sx5 = mcdc.Surface.PlaneX(x=60.0, boundary_condition="vacuum") +sy1 = mcdc.Surface.PlaneY(y=0.0, boundary_condition="reflective") +sy2 = mcdc.Surface.PlaneY(y=10.0) +sy3 = mcdc.Surface.PlaneY(y=50.0) +sy4 = mcdc.Surface.PlaneY(y=60.0) +sy5 = mcdc.Surface.PlaneY(y=100.0, boundary_condition="vacuum") +sz1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective") +sz2 = mcdc.Surface.PlaneZ(z=10.0) +sz3 = mcdc.Surface.PlaneZ(z=30.0) +sz4 = mcdc.Surface.PlaneZ(z=40.0) +sz5 = mcdc.Surface.PlaneZ(z=60.0, boundary_condition="vacuum") + +# Set cells +# Source +mcdc.Cell(region=+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, fill=m) +# Voids +channel_1 = +sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2 +channel_2 = +sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2 +channel_3 = +sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3 +channel_4 = +sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4 +void_channel = channel_1 | channel_2 | channel_3 | channel_4 +mcdc.Cell(region=void_channel, fill=m_void) +# Shield +box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5 +mcdc.Cell(region=box & ~void_channel, fill=m) + +# ====================================================================================== +# Set source +# ====================================================================================== +# The source pulses in t=[0,5] + +mcdc.Source( + x=[0.0, 10.0], + y=[0.0, 10.0], + z=[0.0, 10.0], + isotropic=True, + energy_group=0, +) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +# Tallies +mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100), z=(0.0, 1.0, 60)) +mcdc.TallyMesh(mesh=mesh, scores=["flux"]) + +# Settings +mcdc.settings.N_particle = 1000 +mcdc.settings.N_batch = 2 + +# Techniques +mcdc.simulation.implicit_capture() + +# Run +mcdc.run() diff --git a/examples/kobayashi/process-output.py b/examples/kobayashi/process-output.py new file mode 100644 index 00000000..da2a0f6c --- /dev/null +++ b/examples/kobayashi/process-output.py @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import h5py + +# Load result and grid +with h5py.File("output.h5", "r") as f: + x = f["tallies/mesh_tally_0/grid/x"][:] + y = f["tallies/mesh_tally_0/grid/y"][:] + z = f["tallies/mesh_tally_0/grid/z"][:] + + phi = f["tallies/mesh_tally_0/flux/mean"][:] + phi_sd = f["tallies/mesh_tally_0/flux/sdev"][:] + +# The 2D grid for Z-scan plots +X, Y = np.meshgrid(y, x) + +# Normalization over all Z slices +norm_mean = colors.Normalize(vmin=phi.min(), vmax=phi.max()) +norm_sdev = colors.Normalize(vmin=phi_sd.min(), vmax=phi_sd.max()) + +# Z-scan loop +Nz = len(z) - 1 +for i in range(Nz): + # Get the mean and sdev values for current Z-slice + # The indexing is [x,y,z] + mean = phi[:, :, i] + sdev = phi_sd[:, :, i] + + # Plot mean and sdev side-by-side + fig, ax = plt.subplots(1, 2) + plot_mean = ax[0].pcolormesh(X, Y, mean, norm=norm_mean) + plot_sdev = ax[1].pcolormesh(X, Y, sdev, norm=norm_sdev) + + # Colorbar + fig.colorbar(plot_mean, ax=ax[0]) + fig.colorbar(plot_sdev, ax=ax[1]) + + # Formats + fig.suptitle(f"Flux within Z = [{z[i]}, {z[i+1]}]") + ax[0].set_title("Mean") + ax[0].set_xlabel("$y$ [cm]") + ax[0].set_ylabel("$x$ [cm]") + ax[1].set_title("Std. Dev.") + ax[1].set_xlabel("$y$ [cm]") + ax[1].set_ylabel("$x$ [cm]") + + plt.show() From 8888d4f50d6152b7ba3540fac700d597e04dffb7 Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Tue, 23 Dec 2025 17:49:57 -0800 Subject: [PATCH 20/22] electron related mcdc_get and mcdc_set files --- mcdc/mcdc_get/__init__.py | 10 + .../electron_bremsstrahlung_reaction.py | 3 + .../electron_elastic_scattering_reaction.py | 3 + mcdc/mcdc_get/electron_excitation_reaction.py | 3 + mcdc/mcdc_get/electron_ionization_reaction.py | 32 ++ mcdc/mcdc_get/element.py | 322 ++++++++++++++++++ .../neutron_inelastic_scattering_reaction.py | 12 +- mcdc/mcdc_set/__init__.py | 10 + .../electron_bremsstrahlung_reaction.py | 3 + .../electron_elastic_scattering_reaction.py | 3 + mcdc/mcdc_set/electron_excitation_reaction.py | 3 + mcdc/mcdc_set/electron_ionization_reaction.py | 32 ++ mcdc/mcdc_set/element.py | 322 ++++++++++++++++++ .../neutron_inelastic_scattering_reaction.py | 24 +- .../tabulated_energy_angle_distribution.py | 4 +- 15 files changed, 756 insertions(+), 30 deletions(-) create mode 100644 mcdc/mcdc_get/electron_bremsstrahlung_reaction.py create mode 100644 mcdc/mcdc_get/electron_elastic_scattering_reaction.py create mode 100644 mcdc/mcdc_get/electron_excitation_reaction.py create mode 100644 mcdc/mcdc_get/electron_ionization_reaction.py create mode 100644 mcdc/mcdc_get/element.py create mode 100644 mcdc/mcdc_set/electron_bremsstrahlung_reaction.py create mode 100644 mcdc/mcdc_set/electron_elastic_scattering_reaction.py create mode 100644 mcdc/mcdc_set/electron_excitation_reaction.py create mode 100644 mcdc/mcdc_set/electron_ionization_reaction.py create mode 100644 mcdc/mcdc_set/element.py diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py index 31cde0a6..8e1682bf 100644 --- a/mcdc/mcdc_get/__init__.py +++ b/mcdc/mcdc_get/__init__.py @@ -44,6 +44,16 @@ import mcdc.mcdc_get.tabulated_energy_angle_distribution as tabulated_energy_angle_distribution +import mcdc.mcdc_get.element as element + +import mcdc.mcdc_get.electron_bremsstrahlung_reaction as electron_bremsstrahlung_reaction + +import mcdc.mcdc_get.electron_elastic_scattering_reaction as electron_elastic_scattering_reaction + +import mcdc.mcdc_get.electron_excitation_reaction as electron_excitation_reaction + +import mcdc.mcdc_get.electron_ionization_reaction as electron_ionization_reaction + import mcdc.mcdc_get.gpu_meta as gpu_meta import mcdc.mcdc_get.native_material as native_material diff --git a/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py b/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_elastic_scattering_reaction.py b/mcdc/mcdc_get/electron_elastic_scattering_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_elastic_scattering_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_excitation_reaction.py b/mcdc/mcdc_get/electron_excitation_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_excitation_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_ionization_reaction.py b/mcdc/mcdc_get/electron_ionization_reaction.py new file mode 100644 index 00000000..d05a69d9 --- /dev/null +++ b/mcdc/mcdc_get/electron_ionization_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def binding_energy(index, electron_ionization_reaction, data): + offset = electron_ionization_reaction["binding_energy_offset"] + return data[offset + index] + + +@njit +def binding_energy_all(electron_ionization_reaction, data): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + return data[start:end] + + +@njit +def binding_energy_last(electron_ionization_reaction, data): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + return data[end - 1] + + +@njit +def binding_energy_chunk(start, length, electron_ionization_reaction, data): + start += electron_ionization_reaction["binding_energy_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/element.py b/mcdc/mcdc_get/element.py new file mode 100644 index 00000000..33f86e4f --- /dev/null +++ b/mcdc/mcdc_get/element.py @@ -0,0 +1,322 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs_energy_grid(index, element, data): + offset = element["xs_energy_grid_offset"] + return data[offset + index] + + +@njit +def xs_energy_grid_all(element, data): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + return data[start:end] + + +@njit +def xs_energy_grid_last(element, data): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def xs_energy_grid_chunk(start, length, element, data): + start += element["xs_energy_grid_offset"] + end = start + length + return data[start:end] + + +@njit +def total_xs(index, element, data): + offset = element["total_xs_offset"] + return data[offset + index] + + +@njit +def total_xs_all(element, data): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + return data[start:end] + + +@njit +def total_xs_last(element, data): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def total_xs_chunk(start, length, element, data): + start += element["total_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_xs(index, element, data): + offset = element["ionization_xs_offset"] + return data[offset + index] + + +@njit +def ionization_xs_all(element, data): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + return data[start:end] + + +@njit +def ionization_xs_last(element, data): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def ionization_xs_chunk(start, length, element, data): + start += element["ionization_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def elastic_xs(index, element, data): + offset = element["elastic_xs_offset"] + return data[offset + index] + + +@njit +def elastic_xs_all(element, data): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + return data[start:end] + + +@njit +def elastic_xs_last(element, data): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def elastic_xs_chunk(start, length, element, data): + start += element["elastic_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def excitation_xs(index, element, data): + offset = element["excitation_xs_offset"] + return data[offset + index] + + +@njit +def excitation_xs_all(element, data): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + return data[start:end] + + +@njit +def excitation_xs_last(element, data): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def excitation_xs_chunk(start, length, element, data): + start += element["excitation_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def bremsstrahlung_xs(index, element, data): + offset = element["bremsstrahlung_xs_offset"] + return data[offset + index] + + +@njit +def bremsstrahlung_xs_all(element, data): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + return data[start:end] + + +@njit +def bremsstrahlung_xs_last(element, data): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def bremsstrahlung_xs_chunk(start, length, element, data): + start += element["bremsstrahlung_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_reaction_IDs(index, element, data): + offset = element["ionization_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def ionization_reaction_IDs_all(element, data): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + return data[start:end] + + +@njit +def ionization_reaction_IDs_last(element, data): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + return data[end - 1] + + +@njit +def ionization_reaction_IDs_chunk(start, length, element, data): + start += element["ionization_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def elastic_scattering_reaction_IDs(index, element, data): + offset = element["elastic_scattering_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def elastic_scattering_reaction_IDs_all(element, data): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + return data[start:end] + + +@njit +def elastic_scattering_reaction_IDs_last(element, data): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + return data[end - 1] + + +@njit +def elastic_scattering_reaction_IDs_chunk(start, length, element, data): + start += element["elastic_scattering_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def excitation_reaction_IDs(index, element, data): + offset = element["excitation_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def excitation_reaction_IDs_all(element, data): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + return data[start:end] + + +@njit +def excitation_reaction_IDs_last(element, data): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + return data[end - 1] + + +@njit +def excitation_reaction_IDs_chunk(start, length, element, data): + start += element["excitation_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def bremsstrahlung_reaction_IDs(index, element, data): + offset = element["bremsstrahlung_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def bremsstrahlung_reaction_IDs_all(element, data): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + return data[start:end] + + +@njit +def bremsstrahlung_reaction_IDs_last(element, data): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + return data[end - 1] + + +@njit +def bremsstrahlung_reaction_IDs_chunk(start, length, element, data): + start += element["bremsstrahlung_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_subshell_binding_energy(index, element, data): + offset = element["ionization_subshell_binding_energy_offset"] + return data[offset + index] + + +@njit +def ionization_subshell_binding_energy_all(element, data): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + return data[start:end] + + +@njit +def ionization_subshell_binding_energy_last(element, data): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + return data[end - 1] + + +@njit +def ionization_subshell_binding_energy_chunk(start, length, element, data): + start += element["ionization_subshell_binding_energy_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py b/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py index 1677db2b..5cf2f0dd 100644 --- a/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py +++ b/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py @@ -26,9 +26,7 @@ def spectrum_probability_grid_last(neutron_inelastic_scattering_reaction, data): @njit -def spectrum_probability_grid_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def spectrum_probability_grid_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] end = start + length return data[start:end] @@ -51,9 +49,7 @@ def spectrum_probability(index_1, index_2, neutron_inelastic_scattering_reaction @njit -def spectrum_probability_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def spectrum_probability_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["spectrum_probability_offset"] end = start + length return data[start:end] @@ -82,9 +78,7 @@ def energy_spectrum_IDs_last(neutron_inelastic_scattering_reaction, data): @njit -def energy_spectrum_IDs_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def energy_spectrum_IDs_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] end = start + length return data[start:end] diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py index 4df13dfc..5664c21f 100644 --- a/mcdc/mcdc_set/__init__.py +++ b/mcdc/mcdc_set/__init__.py @@ -44,6 +44,16 @@ import mcdc.mcdc_set.tabulated_energy_angle_distribution as tabulated_energy_angle_distribution +import mcdc.mcdc_set.element as element + +import mcdc.mcdc_set.electron_bremsstrahlung_reaction as electron_bremsstrahlung_reaction + +import mcdc.mcdc_set.electron_elastic_scattering_reaction as electron_elastic_scattering_reaction + +import mcdc.mcdc_set.electron_excitation_reaction as electron_excitation_reaction + +import mcdc.mcdc_set.electron_ionization_reaction as electron_ionization_reaction + import mcdc.mcdc_set.gpu_meta as gpu_meta import mcdc.mcdc_set.native_material as native_material diff --git a/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py b/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_elastic_scattering_reaction.py b/mcdc/mcdc_set/electron_elastic_scattering_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_elastic_scattering_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_excitation_reaction.py b/mcdc/mcdc_set/electron_excitation_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_excitation_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_ionization_reaction.py b/mcdc/mcdc_set/electron_ionization_reaction.py new file mode 100644 index 00000000..df76b1b9 --- /dev/null +++ b/mcdc/mcdc_set/electron_ionization_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def binding_energy(index, electron_ionization_reaction, data, value): + offset = electron_ionization_reaction["binding_energy_offset"] + data[offset + index] = value + + +@njit +def binding_energy_all(electron_ionization_reaction, data, value): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + data[start:end] = value + + +@njit +def binding_energy_last(electron_ionization_reaction, data, value): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + data[end - 1] = value + + +@njit +def binding_energy_chunk(start, length, electron_ionization_reaction, data, value): + start += electron_ionization_reaction["binding_energy_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/element.py b/mcdc/mcdc_set/element.py new file mode 100644 index 00000000..11a248a9 --- /dev/null +++ b/mcdc/mcdc_set/element.py @@ -0,0 +1,322 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs_energy_grid(index, element, data, value): + offset = element["xs_energy_grid_offset"] + data[offset + index] = value + + +@njit +def xs_energy_grid_all(element, data, value): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def xs_energy_grid_last(element, data, value): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def xs_energy_grid_chunk(start, length, element, data, value): + start += element["xs_energy_grid_offset"] + end = start + length + data[start:end] = value + + +@njit +def total_xs(index, element, data, value): + offset = element["total_xs_offset"] + data[offset + index] = value + + +@njit +def total_xs_all(element, data, value): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def total_xs_last(element, data, value): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def total_xs_chunk(start, length, element, data, value): + start += element["total_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_xs(index, element, data, value): + offset = element["ionization_xs_offset"] + data[offset + index] = value + + +@njit +def ionization_xs_all(element, data, value): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def ionization_xs_last(element, data, value): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_xs_chunk(start, length, element, data, value): + start += element["ionization_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def elastic_xs(index, element, data, value): + offset = element["elastic_xs_offset"] + data[offset + index] = value + + +@njit +def elastic_xs_all(element, data, value): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def elastic_xs_last(element, data, value): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def elastic_xs_chunk(start, length, element, data, value): + start += element["elastic_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def excitation_xs(index, element, data, value): + offset = element["excitation_xs_offset"] + data[offset + index] = value + + +@njit +def excitation_xs_all(element, data, value): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def excitation_xs_last(element, data, value): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def excitation_xs_chunk(start, length, element, data, value): + start += element["excitation_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def bremsstrahlung_xs(index, element, data, value): + offset = element["bremsstrahlung_xs_offset"] + data[offset + index] = value + + +@njit +def bremsstrahlung_xs_all(element, data, value): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def bremsstrahlung_xs_last(element, data, value): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def bremsstrahlung_xs_chunk(start, length, element, data, value): + start += element["bremsstrahlung_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_reaction_IDs(index, element, data, value): + offset = element["ionization_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def ionization_reaction_IDs_all(element, data, value): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + data[start:end] = value + + +@njit +def ionization_reaction_IDs_last(element, data, value): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_reaction_IDs_chunk(start, length, element, data, value): + start += element["ionization_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def elastic_scattering_reaction_IDs(index, element, data, value): + offset = element["elastic_scattering_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def elastic_scattering_reaction_IDs_all(element, data, value): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + data[start:end] = value + + +@njit +def elastic_scattering_reaction_IDs_last(element, data, value): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def elastic_scattering_reaction_IDs_chunk(start, length, element, data, value): + start += element["elastic_scattering_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def excitation_reaction_IDs(index, element, data, value): + offset = element["excitation_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def excitation_reaction_IDs_all(element, data, value): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + data[start:end] = value + + +@njit +def excitation_reaction_IDs_last(element, data, value): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def excitation_reaction_IDs_chunk(start, length, element, data, value): + start += element["excitation_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def bremsstrahlung_reaction_IDs(index, element, data, value): + offset = element["bremsstrahlung_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def bremsstrahlung_reaction_IDs_all(element, data, value): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + data[start:end] = value + + +@njit +def bremsstrahlung_reaction_IDs_last(element, data, value): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def bremsstrahlung_reaction_IDs_chunk(start, length, element, data, value): + start += element["bremsstrahlung_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_subshell_binding_energy(index, element, data, value): + offset = element["ionization_subshell_binding_energy_offset"] + data[offset + index] = value + + +@njit +def ionization_subshell_binding_energy_all(element, data, value): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + data[start:end] = value + + +@njit +def ionization_subshell_binding_energy_last(element, data, value): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_subshell_binding_energy_chunk(start, length, element, data, value): + start += element["ionization_subshell_binding_energy_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py b/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py index 65509585..81468be2 100644 --- a/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py +++ b/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py @@ -4,9 +4,7 @@ @njit -def spectrum_probability_grid( - index, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_grid(index, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] data[offset + index] = value @@ -28,18 +26,14 @@ def spectrum_probability_grid_last(neutron_inelastic_scattering_reaction, data, @njit -def spectrum_probability_grid_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_grid_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] end = start + length data[start:end] = value @njit -def spectrum_probability_vector( - index_1, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_vector(index_1, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_offset"] stride = neutron_inelastic_scattering_reaction["N_spectrum"] start = offset + index_1 * stride @@ -48,18 +42,14 @@ def spectrum_probability_vector( @njit -def spectrum_probability( - index_1, index_2, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability(index_1, index_2, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_offset"] stride = neutron_inelastic_scattering_reaction["N_spectrum"] data[offset + index_1 * stride + index_2] = value @njit -def spectrum_probability_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["spectrum_probability_offset"] end = start + length data[start:end] = value @@ -88,9 +78,7 @@ def energy_spectrum_IDs_last(neutron_inelastic_scattering_reaction, data, value) @njit -def energy_spectrum_IDs_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def energy_spectrum_IDs_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] end = start + length data[start:end] = value diff --git a/mcdc/mcdc_set/tabulated_energy_angle_distribution.py b/mcdc/mcdc_set/tabulated_energy_angle_distribution.py index 97040280..e25747d1 100644 --- a/mcdc/mcdc_set/tabulated_energy_angle_distribution.py +++ b/mcdc/mcdc_set/tabulated_energy_angle_distribution.py @@ -171,9 +171,7 @@ def cosine_offset__last(tabulated_energy_angle_distribution, data, value): @njit -def cosine_offset__chunk( - start, length, tabulated_energy_angle_distribution, data, value -): +def cosine_offset__chunk(start, length, tabulated_energy_angle_distribution, data, value): start += tabulated_energy_angle_distribution["cosine_offset__offset"] end = start + length data[start:end] = value From 155ee400594d09267aaeb124c2ffa52eb121f54d Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Tue, 23 Dec 2025 17:52:36 -0800 Subject: [PATCH 21/22] some minor fixes before proceeding --- mcdc/object_/material.py | 1 + mcdc/object_/reaction.py | 33 ++++++++++++++++++--------------- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/mcdc/object_/material.py b/mcdc/object_/material.py index cfbd7d47..a608b719 100644 --- a/mcdc/object_/material.py +++ b/mcdc/object_/material.py @@ -11,6 +11,7 @@ from mcdc.constant import MATERIAL, MATERIAL_MG from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.nuclide import Nuclide +from mcdc.object_.element import Element from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index 2a9684c2..529a3fd0 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -23,6 +23,7 @@ REACTION_ELECTRON_IONIZATION, REFERENCE_FRAME_COM, REFERENCE_FRAME_LAB, + MU_CUTOFF, ) from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.distribution import ( @@ -50,7 +51,7 @@ class ReactionBase(ObjectPolymorphic): # MT: int xs: NDArray[float64] - xs_offset_: int # "xs_offset" ir reserved for "xs" + xs_offset_: int # "xs_offset" is reserved for "xs" reference_frame: int def __init__(self, type_, MT, xs, xs_offset, reference_frame): @@ -305,8 +306,9 @@ class ReactionElectronIonization(ReactionBase): label: str = "electron_ionization_reaction" # N_subshell: int - subshell_xs: list[DataBase] - subshell_product: list[DistributionBase] + subshell_xs: DistributionBase + binding_energy: NDArray[float64] + subshell_product: DistributionBase def __init__( self, @@ -315,6 +317,7 @@ def __init__( xs_offset, reference_frame, subshell_xs, + binding_energy, subshell_product, ): type_ = REACTION_ELECTRON_IONIZATION @@ -322,6 +325,7 @@ def __init__( self.N_subshell = len(subshell_xs) self.subshell_xs = subshell_xs + self.binding_energy = binding_energy self.subshell_product = subshell_product @classmethod @@ -333,6 +337,7 @@ def from_h5_group(cls, h5_group): subshell_xs = [] subshell_product = [] + binding_energy = [] for name in subshell_names: subshell = subshells[name] @@ -343,6 +348,7 @@ def from_h5_group(cls, h5_group): subshell["xs"][()] ) subshell_xs.append(xs_sub) + binding_energy.append(subshell[0].x[0]) # Secondary electron energy distribution product = subshell["product"] @@ -360,6 +366,7 @@ def from_h5_group(cls, h5_group): xs_offset, reference_frame, subshell_xs, + binding_energy, subshell_product, ) @@ -385,9 +392,8 @@ def __repr__(self): class ReactionElectronElasticScattering(ReactionBase): label: str = "electron_elastic_scattering_reaction" # - mu_cut: float xs_large: DataBase - mu_table: DistributionMultiTable + mu: DistributionMultiTable def __init__( self, @@ -395,15 +401,14 @@ def __init__( xs, xs_offset, reference_frame, - mu_cut, xs_large, - mu_table, + mu, ): type_ = REACTION_ELECTRON_ELASTIC_SCATTERING super().__init__(type_, MT, xs, xs_offset, reference_frame) - self.mu_cut = mu_cut + self.mu_cut = MU_CUTOFF self.xs_large = xs_large - self.mu_table = mu_table + self.mu = mu @classmethod def from_h5_group(cls, h5_group): @@ -411,7 +416,6 @@ def from_h5_group(cls, h5_group): # Large angle data large_angle = h5_group["large_angle"] - mu_cut = float(large_angle.attrs.get("mu_cut", 0.999999)) xs_large = DataTable( large_angle["xs_energy"][()], @@ -419,7 +423,7 @@ def from_h5_group(cls, h5_group): ) mu_g = large_angle["scattering_cosine"] - mu_table = DistributionMultiTable( + mu = DistributionMultiTable( mu_g["energy_grid"][()], mu_g["energy_offset"][()], mu_g["value"][()], @@ -431,9 +435,8 @@ def from_h5_group(cls, h5_group): xs, xs_offset, reference_frame, - mu_cut, xs_large, - mu_table, + mu, ) def __repr__(self): @@ -441,8 +444,8 @@ def __repr__(self): text += f" - Mu cut: {self.mu_cut}\n" text += f" - Large angle XS: DataTable [ID: {self.xs_large.ID}]\n" text += ( - f" - Scattering cosine: {distribution.decode_type(self.mu_table.type)} " - f"[ID: {self.mu_table.ID}]\n" + f" - Scattering cosine: {distribution.decode_type(self.mu.type)} " + f"[ID: {self.mu.ID}]\n" ) return text From 69e7b16d4b42933bb0641de551424cc2d7a62081 Mon Sep 17 00:00:00 2001 From: Melek Derman <48313913+melekderman@users.noreply.github.com> Date: Wed, 24 Dec 2025 02:50:03 -0800 Subject: [PATCH 22/22] add element support to material.py --- mcdc/object_/material.py | 183 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 182 insertions(+), 1 deletion(-) diff --git a/mcdc/object_/material.py b/mcdc/object_/material.py index a608b719..6d9a5acb 100644 --- a/mcdc/object_/material.py +++ b/mcdc/object_/material.py @@ -14,6 +14,7 @@ from mcdc.object_.element import Element from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error +from mcdc.object_.util import ISOTOPIC_ABUNDANCE # ====================================================================================== # Material base class @@ -63,17 +64,22 @@ class Material(MaterialBase): # Annotations for Numba mode label: str = "native_material" non_numba: list[str] = ["nuclide_composition"] + non_numba_element: list[str] = ["element_composition"] # nuclide_composition: dict[Nuclide, float] + element_composition: dict[Element, float] # nuclides: list[Nuclide] + elements: list[Element] nuclide_densities: NDArray[float64] + element_densities: NDArray[float64] def __init__( self, name: str = "", nuclide_composition: dict[str, float] = {}, - temperature: float = 293.6, + element_composition: dict[str, float] = {}, + temperature: float = 293.6 ): type_ = MATERIAL super().__init__(type_, name) @@ -84,15 +90,73 @@ def __init__( # Dictionary connecting nuclides to respective densities self.nuclide_composition = {} + # Dictionary connecting elements to respective densities + self.element_composition = {} + # Numba representation of nuclide_composition self.nuclides = [] self.nuclide_densities = np.zeros(len(nuclide_composition)) + # Numba representation of element_composition + self.elements = [] + self.element_densities = np.zeros(len(element_composition)) + + # ============================================================================== + # Argument checks + # ============================================================================== + + # Check that either nuclide or element composition is provided, not both + if len(nuclide_composition) > 0 and len(element_composition) > 0: + print_error("Cannot specify both nuclide_composition and element_composition") + + if len(nuclide_composition) == 0 and len(element_composition) == 0: + print_error("Must specify either nuclide_composition or element_composition") + # Check if library directory is set lib_dir = os.getenv("MCDC_LIB") if lib_dir is None: print_error("Environment variable MCDC_LIB is not set") + # ============================================================================== + # Element-based construction + # ============================================================================== + + # Loop over the items in the composition + for i, (key, value) in enumerate(element_composition.items()): + element_name = key + element_density = value + temperature = 0.0 # Elements are always at 0 K + + # Check if element is available in the library + file_name = f"{element_name}.h5" + if not file_name in os.listdir(lib_dir): + print_error( + f"Element {element_name} is not available in the library" + ) + + # Check if element is already created + found = False + for element in simulation.elements: + if ( + element.name == element_name + ): + found = True + break + + # Create the element to objects if needed + if not found: + element = Element(element_name) + + # Register the element composition + self.elements.append(element) + self.element_densities[i] = element_density + self.element_composition[element] = element_density + + + # ============================================================================== + # Nuclide-based construction + # ============================================================================== + # Loop over the items in the composition for i, (key, value) in enumerate(nuclide_composition.items()): nuclide_name = key @@ -131,17 +195,27 @@ def __init__( if nuclide.fissionable: self.fissionable = True + # Create undefined compositions + if self.nuclides == []: + set_nuclides_from_elements(self) + if self.elements == []: + set_elements_from_nuclides(self) + def __repr__(self): text = super().__repr__() text += f" - Temperature: {self.temperature} K\n" text += f" - Nuclide composition [atoms/barn-cm]\n" for nuclide in self.nuclide_composition.keys(): text += f" - {nuclide.name:<5} | {self.nuclide_composition[nuclide]}\n" + text += f" - Element composition [atoms/barn-cm]\n" + for element in self.element_composition.keys(): + text += f" - {element.name:<5} | {self.element_composition[element]}\n" return text # Currently supported temperatures TEMPERATURES = [0.1, 233.15, 273.15, 293.6, 600.0, 900.0, 1200.0, 2500.0] +TEMPERATURES_ELEMENT = [0.0] # ====================================================================================== @@ -322,3 +396,110 @@ def __repr__(self): text += f" - speed {print_1d_array(self.mgxs_speed)}\n" text += f" - lambda {print_1d_array(self.mgxs_decay_rate)}\n" return text + + +# ====================================================================================== +# Helpers +# ====================================================================================== + +def set_nuclides_from_elements(material): + material.nuclides = [] + material.nuclide_densities = [] + material.nuclide_composition = {} + + # Check if library directory is set + lib_dir = os.getenv("MCDC_LIB") + if lib_dir is None: + print_error("Environment variable MCDC_LIB is not set") + + # Get supported temperature + nearest_temperature = min( + TEMPERATURES, key=lambda x: abs(x - material.temperature) + ) + + for element, element_density in material.element_composition.items(): + # To make sure that the abundance is normalized + norm = 0.0 + for (_, abundance) in ISOTOPIC_ABUNDANCE[element.name].items(): + norm += abundance + + # Loop over the nuclide composition + for nuclide_name, abundance in ISOTOPIC_ABUNDANCE[element.name].items(): + # Check if nuclide-temperature is available in the library + file_name = f"{nuclide_name}-{nearest_temperature}K.h5" + if not file_name in os.listdir(lib_dir): + print_error( + f"Nuclide {nuclide_name} at temperature {nearest_temperature} K is not available in the library" + ) + + # Check if nuclide is already created + found = False + for nuclide in simulation.nuclides: + if ( + nuclide.name == nuclide_name + and nearest_temperature == nuclide.temperature + ): + found = True + break + + # Create the nuclide to objects if needed + if not found: + nuclide = Nuclide(nuclide_name, nearest_temperature) + + # Calculate nuclide density + nuclide_density = element_density * abundance / norm + + # Register the nuclide composition + material.nuclides.append(nuclide) + material.nuclide_densities.append(nuclide_density) + material.nuclide_composition[nuclide] = nuclide_density + + # Some flags + if nuclide.fissionable: + material.fissionable = True + + material.nuclide_densities = np.array(material.nuclide_densities) + + +def set_elements_from_nuclides(material): + material.elements = [] + material.element_densities = [] + material.element_composition = {} + + # Get the list of the element names + element_names = [] + for nuclide in material.nuclides: + name = nuclide.name[:2] + if name[1].isdigit(): + name = name[0] + if name not in element_names: + element_names.append(name) + material.element_densities = np.zeros(len(element_names)) + + # Iterate over all named elements + for i, element_name in enumerate(element_names): + # Check if element is already created + found = False + for element in simulation.elements: + if element.name == element_name: + found = True + break + + # Create the element to objects if needed + if not found: + element = Element(element_name) + + material.elements.append(element) + + # Iterate over all nuclides to get the total density + density = 0.0 + for nuclide, nuclide_density in material.nuclide_composition.items(): + # Skip if non-isotope + if nuclide.name[: len(element_name)] != element_name: + continue + + # Accumulate density + density += material.nuclide_composition[nuclide] + + material.element_densities[i] = density + material.element_composition[element] = density