diff --git a/tools/data_library_generator/electron/README.md b/tools/data_library_generator/electron/README.md new file mode 100644 index 000000000..7c6724207 --- /dev/null +++ b/tools/data_library_generator/electron/README.md @@ -0,0 +1,90 @@ +# MC/DC Electron Data Library Generator +Converts EPRDATA14 ACE-format electron/photon/relaxation data into MC/DC's +per-element HDF5 format for continuous-energy electron transport. + +## Prerequisites + +- Installing ACEtk from source: [link](https://github.com/njoy/ACEtk) +- Dependencies: `pip install h5py numpy tqdm` +- You need the EPRDATA14 library (available from [LANL Nuclear Data](https://nucleardata.lanl.gov/ace/eprdata14)). + +## Environment Variables +| Variable | Description | +|------------------------|----------------------------------------------------| +| `MCDC_ACELIB_ELECTRON` | Path to the EPRDATA14 data file. | +| `MCDC_LIB_ELECTRON` | Path to the output directory for MC/DC HDF5 files. | + +## Usage +```bash +export MCDC_ACELIB_ELECTRON=/path/to/eprdata14/eprdata14/eprdata14 +export MCDC_LIB_ELECTRON=/path/to/mcdc/electron/library + +python generate.py # Convert only missing elements +python generate.py --rewrite # Regenerate all files +python generate.py --verbose # Print detailed per-element info +``` + +## What it Does +For each element (Z=1 to Z=100) in the EPRDATA14 library, the generator: +1. Loads all elemental tables from the single concatenated EPRDATA14 file. +2. Extracts the principal cross section energy grid and pointwise cross sections + (elastic, bremsstrahlung, excitation, ionization) for all reaction channels. +3. Extracts tabulated elastic angular distributions (cosine CDFs) per incident energy (MT-528). +4. Extracts excitation energy loss as a function of incident energy (MT-527). +5. Extracts bremsstrahlung outgoing photon energy distributions per incident energy (MT-526). +6. Extracts per-subshell ionization cross sections and knock-on electron + energy distributions (MT-534 for K-shell, MT-535+ for higher shells). +7. Extracts atomic relaxation (fluorescence and Auger) transition data per subshell. +8. Writes a single HDF5 file per element (e.g., `Al.h5`). + +## Output HDF5 Schema +``` +.h5 +├── element_symbol (string) +├── atomic_number (int) +├── atomic_weight_ratio (float) +├── electron_reactions/ +│ ├── xs_energy_grid (1-D array, MeV) +│ ├── elastic_scattering/ +│ │ └── MT-528/ +│ │ ├── xs (1-D array, barns) +│ │ └── angular_cosine_distribution/ +│ │ ├── energy (1-D array, MeV) +│ │ ├── offset (1-D array, int) +│ │ ├── cosine (1-D array) +│ │ └── cdf (1-D array) +│ ├── excitation/ +│ │ └── MT-527/ +│ │ ├── xs (1-D array, barns) +│ │ └── energy_loss/ +│ │ ├── energy (1-D array, MeV) +│ │ └── excitation_energy_loss (1-D array, MeV) +│ ├── bremsstrahlung/ +│ │ └── MT-526/ +│ │ ├── xs (1-D array, barns) +│ │ └── energy_distribution/ +│ │ ├── energy (1-D array, MeV) +│ │ ├── offset (1-D array, int) +│ │ ├── energy_out (1-D array, MeV) +│ │ └── cdf (1-D array) +│ └── ionization/ +│ └── MT-534/ (K-shell; MT-535, MT-536, ... for higher shells) +│ ├── xs (1-D array, barns) +│ ├── binding_energy (float, MeV) +│ └── energy_distribution/ +│ ├── energy (1-D array, MeV) +│ ├── offset (1-D array, int) +│ ├── energy_out (1-D array, MeV) +│ └── cdf (1-D array) +└── atomic_relaxation/ + └── MT-534/ (K-shell; MT-535, MT-536, ... for higher shells) + ├── number_of_transitions (int) + ├── primary_designator (1-D array, int) + ├── secondary_designator (1-D array, int) + ├── energy (1-D array, MeV) + └── probability (1-D array) +``` + +## See Also +- [Continuous Energy Theory Guide](../../docs/source/theory/cont_energy.rst) +- [Installation Guide — CE Library Configuration](../../docs/source/install.rst) \ No newline at end of file diff --git a/tools/data_library_generator/electron/generate.py b/tools/data_library_generator/electron/generate.py new file mode 100644 index 000000000..6de516b96 --- /dev/null +++ b/tools/data_library_generator/electron/generate.py @@ -0,0 +1,247 @@ +import ACEtk +import argparse +import h5py +import numpy as np +import os + +from tqdm import tqdm + +#### + +import util +from util import print_error, print_note + +parser = argparse.ArgumentParser(description="MC/DC electron data generator") +parser.add_argument("--rewrite", dest="rewrite", action="store_true", default=False) +parser.add_argument("--verbose", dest="verbose", action="store_true", default=False) +args, unargs = parser.parse_known_args() +rewrite = args.rewrite +verbose = args.verbose + +# Directories +output_dir = os.getenv("MCDC_LIB_ELECTRON") +ace_file = os.getenv("MCDC_ACELIB_ELECTRON") +if output_dir is None: + print_error("Environment variable $MCDC_LIB_ELECTRON is not set") +if ace_file is None: + print_error("Environment variable $MCDC_ACELIB_ELECTRON is not set") +# Create output directory if needed +os.makedirs(output_dir, exist_ok=True) +print(f"\nACE file : {ace_file}") +print(f"Output dir : {output_dir}\n") + +# Load all tables from the concatenated EPRDATA14 file +print("Loading EPRDATA14 tables...") +all_tables = ACEtk.PhotoatomicTable.from_concatenated_file(ace_file) +table_map = {t.zaid: t for t in all_tables} +print(f"Loaded {len(table_map)} tables\n") + +# Select target entries +target_entries = [] +for zaid in table_map: + if not zaid.endswith(".14p"): + continue + + Z = util.decode_epr_zaid(zaid) + symbol = util.Z_TO_SYMBOL[Z] + mcdc_name = f"{symbol}.h5" + + if not rewrite and os.path.exists(f"{output_dir}/{mcdc_name}"): + continue + + target_entries.append((zaid, Z, symbol, mcdc_name)) + +# Loop over all elements +pbar = tqdm( + target_entries, + disable=verbose, + bar_format="{l_bar}{bar}| {n_fmt}/{total_fmt}{postfix}", +) +for zaid, Z, symbol, mcdc_name in pbar: + + if not rewrite and os.path.exists(f"{output_dir}/{mcdc_name}"): + continue + + if verbose: + print("\n" + "=" * 80 + "\n") + print(f"Create {mcdc_name} from {zaid}\n") + pbar.set_postfix_str(f"{mcdc_name}") + + # Load ACE table + ace_table = table_map[zaid] + + # Create MC/DC file + file = h5py.File(f"{output_dir}/{mcdc_name}", "w") + + # ================================================================================== + # Basic properties + # ================================================================================== + + header = ace_table.header + file.attrs["source_title"] = header.title + file.attrs["source_date"] = header.date + + file.create_dataset("element_symbol", data=symbol) + file.create_dataset("atomic_number", data=Z) + file.create_dataset("atomic_weight_ratio", data=ace_table.atomic_weight_ratio) + + # ================================================================================== + # Reaction groups + # ================================================================================== + # Elastic scattering : MT-528 + # Excitation : MT-527 + # Bremsstrahlung : MT-526 + # Ionization : MT-534 (K), MT-535 (L1), MT-536 (L2), ... + + reactions = file.create_group("electron_reactions") + + elastic_group = reactions.create_group("elastic_scattering") + excitation_group = reactions.create_group("excitation") + bremsstrahlung_group = reactions.create_group("bremsstrahlung") + ionization_group = reactions.create_group("ionization") + + elastic_group.create_group("MT-528").attrs["MT"] = 528 + excitation_group.create_group("MT-527").attrs["MT"] = 527 + bremsstrahlung_group.create_group("MT-526").attrs["MT"] = 526 + + subsh_block = ace_table.electron_subshell_block + N_subshells = subsh_block.number_electron_subshells + ionization_MTs = [534 + i for i in range(N_subshells)] + + for MT in ionization_MTs: + ionization_group.create_group(f"MT-{MT:03}").attrs["MT"] = MT + + if verbose: + print(f" Reaction group MTs") + print(f" - Elastic scattering MT : [528]") + print(f" - Excitation MT : [527]") + print(f" - Bremsstrahlung MT : [526]") + print(f" - Ionization MTs : {ionization_MTs}") + + # ================================================================================== + # Cross sections + # ================================================================================== + + xs0_block = ace_table.electron_cross_section_block + + xs_energy = np.array(xs0_block.energies) + dataset = reactions.create_dataset("xs_energy_grid", data=xs_energy) + dataset.attrs["unit"] = "MeV" + + xs = elastic_group.create_dataset("MT-528/xs", data=np.array(xs0_block.elastic)) + xs.attrs["unit"] = "barns" + + xs = excitation_group.create_dataset( + "MT-527/xs", data=np.array(xs0_block.excitation) + ) + xs.attrs["unit"] = "barns" + + xs = bremsstrahlung_group.create_dataset( + "MT-526/xs", data=np.array(xs0_block.bremsstrahlung) + ) + xs.attrs["unit"] = "barns" + + for i, MT in enumerate(ionization_MTs): + xs = ionization_group.create_dataset( + f"MT-{MT:03}/xs", data=np.array(xs0_block.electroionisation(i + 1)) + ) + xs.attrs["unit"] = "barns" + + # ================================================================================== + # Elastic angular distribution + # ================================================================================== + + angle_group = elastic_group.create_group("MT-528/angular_cosine_distribution") + util.load_elastic_angular_distribution( + ace_table.electron_elastic_angular_distribution_block, angle_group + ) + + # ================================================================================== + # Excitation energy loss + # ================================================================================== + + excit_block = ace_table.electron_excitation_energy_loss_block + excit_group = excitation_group.create_group("MT-527/energy_loss") + + dataset = excit_group.create_dataset("energy", data=np.array(excit_block.energies)) + dataset.attrs["unit"] = "MeV" + + dataset = excit_group.create_dataset( + "excitation_energy_loss", data=np.array(excit_block.excitation_energy_loss) + ) + dataset.attrs["unit"] = "MeV" + + # ================================================================================== + # Bremsstrahlung energy distribution + # ================================================================================== + + brems_group = bremsstrahlung_group.create_group("MT-526/energy_distribution") + util.load_bremsstrahlung( + ace_table.bremsstrahlung_energy_distribution_block, brems_group + ) + + # ================================================================================== + # Ionization: binding energy and knock-on electron energy distributions + # ================================================================================== + + for i, MT in enumerate(ionization_MTs): + idx = i + 1 + MT_group = ionization_group[f"MT-{MT:03}"] + + dataset = MT_group.create_dataset( + "binding_energy", data=subsh_block.binding_energy(idx) + ) + dataset.attrs["unit"] = "MeV" + + energy_dist_group = MT_group.create_group("energy_distribution") + util.load_electroionization_subshell( + ace_table.electroionisation_energy_distribution_block(idx), + energy_dist_group, + ) + + # ================================================================================== + # Atomic relaxation (subshell transition data) + # ================================================================================== + + relax_block = ace_table.subshell_transition_data_block + relaxation_group = file.create_group("atomic_relaxation") + + for i, MT in enumerate(ionization_MTs): + idx = i + 1 + td = relax_block.transition_data(idx) + MT_group = relaxation_group.create_group(f"MT-{MT:03}") + MT_group.attrs["MT"] = MT + + N_transitions = td.number_transitions + MT_group.create_dataset("number_of_transitions", data=N_transitions) + + if N_transitions > 0: + primary_designators = [] + secondary_designators = [] + energies = [] + probabilities = [] + for j in range(N_transitions): + jdx = j + 1 + t = td.transition(jdx) + primary_designators.append(td.primary_designator(jdx)) + secondary_designators.append(td.secondary_designator(jdx)) + energies.append(td.energy(jdx)) + probabilities.append(td.probability(jdx)) + + MT_group.create_dataset( + "primary_designator", data=np.array(primary_designators) + ) + MT_group.create_dataset( + "secondary_designator", data=np.array(secondary_designators) + ) + dataset = MT_group.create_dataset("energy", data=np.array(energies)) + dataset.attrs["unit"] = "MeV" + MT_group.create_dataset("probability", data=np.array(probabilities)) + + # ================================================================================== + # Finalize + # ================================================================================== + + file.close() + +print("") diff --git a/tools/data_library_generator/electron/util.py b/tools/data_library_generator/electron/util.py new file mode 100644 index 000000000..91de4264c --- /dev/null +++ b/tools/data_library_generator/electron/util.py @@ -0,0 +1,207 @@ +import ACEtk +import h5py +import numpy as np + + +def print_error(message): + print(f"\n [ERROR]: {message}\n") + exit() + + +def print_note(message): + print(f"\n [NOTE]: {message}\n") + + +def decode_epr_zaid(name: str): + """ + Decode an EPR ACE ZAID like '1000.14p' into atomic number Z. + All EPR tables are elemental: ZAID = Z * 1000 (A = 000). + Returns Z. + """ + zaid_str, _ = name.split(".") + Z = int(zaid_str) // 1000 + return Z + + +def load_elastic_angular_distribution(block, h5_group: h5py.Group): + """ + Load elastic angular distribution block into HDF5 group. + Stores incident energy grid, cosine grid, and CDF per energy. + Note: ACEtk provides CDF only (no PDF) for electron elastic angular data. + """ + energies = np.array(block.energies) + dataset = h5_group.create_dataset("energy", data=energies) + dataset.attrs["unit"] = "MeV" + + NE = len(energies) + offset = np.zeros(NE, dtype=int) + cosines = [] + cdf = [] + for i, dist in enumerate(block.distributions): + offset[i] = len(cosines) + cosines.extend(dist.cosines) + cdf.extend(dist.cdf) + + h5_group.create_dataset("offset", data=offset) + h5_group.create_dataset("cosine", data=np.array(cosines)) + h5_group.create_dataset("cdf", data=np.array(cdf)) + + +def load_bremsstrahlung(block, h5_group: h5py.Group): + """ + Load bremsstrahlung energy distribution block into HDF5 group. + Stores outgoing photon energy CDF for each incident electron energy. + Note: ACEtk provides CDF only (no PDF) for bremsstrahlung distributions. + """ + energies = np.array(block.energies) + dataset = h5_group.create_dataset("energy", data=energies) + dataset.attrs["unit"] = "MeV" + + NE = len(energies) + offset = np.zeros(NE, dtype=int) + energy_out = [] + cdf = [] + for i, dist in enumerate(block.distributions): + offset[i] = len(energy_out) + energy_out.extend(dist.outgoing_energies) + cdf.extend(dist.cdf) + + h5_group.create_dataset("offset", data=offset) + dataset = h5_group.create_dataset("energy_out", data=np.array(energy_out)) + dataset.attrs["unit"] = "MeV" + h5_group.create_dataset("cdf", data=np.array(cdf)) + + +def load_electroionization_subshell(block, h5_group: h5py.Group): + """ + Load electroionization energy distribution for a single subshell. + Stores outgoing (knock-on) electron energy CDF per incident energy. + Note: ACEtk provides CDF only (no PDF) for electroionization distributions. + """ + energies = np.array(block.energies) + dataset = h5_group.create_dataset("energy", data=energies) + dataset.attrs["unit"] = "MeV" + + NE = len(energies) + offset = np.zeros(NE, dtype=int) + energy_out = [] + cdf = [] + for i, dist in enumerate(block.distributions): + offset[i] = len(energy_out) + energy_out.extend(dist.outgoing_energies) + cdf.extend(dist.cdf) + + h5_group.create_dataset("offset", data=offset) + dataset = h5_group.create_dataset("energy_out", data=np.array(energy_out)) + dataset.attrs["unit"] = "MeV" + h5_group.create_dataset("cdf", data=np.array(cdf)) + + +# ============================================================================= +# Constants +# ============================================================================= + +SYMBOL_TO_Z = { + "H": 1, + "He": 2, + "Li": 3, + "Be": 4, + "B": 5, + "C": 6, + "N": 7, + "O": 8, + "F": 9, + "Ne": 10, + "Na": 11, + "Mg": 12, + "Al": 13, + "Si": 14, + "P": 15, + "S": 16, + "Cl": 17, + "Ar": 18, + "K": 19, + "Ca": 20, + "Sc": 21, + "Ti": 22, + "V": 23, + "Cr": 24, + "Mn": 25, + "Fe": 26, + "Co": 27, + "Ni": 28, + "Cu": 29, + "Zn": 30, + "Ga": 31, + "Ge": 32, + "As": 33, + "Se": 34, + "Br": 35, + "Kr": 36, + "Rb": 37, + "Sr": 38, + "Y": 39, + "Zr": 40, + "Nb": 41, + "Mo": 42, + "Tc": 43, + "Ru": 44, + "Rh": 45, + "Pd": 46, + "Ag": 47, + "Cd": 48, + "In": 49, + "Sn": 50, + "Sb": 51, + "Te": 52, + "I": 53, + "Xe": 54, + "Cs": 55, + "Ba": 56, + "La": 57, + "Ce": 58, + "Pr": 59, + "Nd": 60, + "Pm": 61, + "Sm": 62, + "Eu": 63, + "Gd": 64, + "Tb": 65, + "Dy": 66, + "Ho": 67, + "Er": 68, + "Tm": 69, + "Yb": 70, + "Lu": 71, + "Hf": 72, + "Ta": 73, + "W": 74, + "Re": 75, + "Os": 76, + "Ir": 77, + "Pt": 78, + "Au": 79, + "Hg": 80, + "Tl": 81, + "Pb": 82, + "Bi": 83, + "Po": 84, + "At": 85, + "Rn": 86, + "Fr": 87, + "Ra": 88, + "Ac": 89, + "Th": 90, + "Pa": 91, + "U": 92, + "Np": 93, + "Pu": 94, + "Am": 95, + "Cm": 96, + "Bk": 97, + "Cf": 98, + "Es": 99, + "Fm": 100, +} + +Z_TO_SYMBOL = {value: key for key, value in SYMBOL_TO_Z.items()} diff --git a/tools/data_library_generator/README.md b/tools/data_library_generator/neutron/README.md similarity index 81% rename from tools/data_library_generator/README.md rename to tools/data_library_generator/neutron/README.md index 43065cd76..91ac4b18e 100644 --- a/tools/data_library_generator/README.md +++ b/tools/data_library_generator/neutron/README.md @@ -5,24 +5,22 @@ for continuous-energy neutron transport. ## Prerequisites -```bash -pip install ACEtk h5py numpy tqdm -``` - -You need a collection of ACE files (e.g., from NJOY or an ENDF/B distribution). +- Installing ACEtk from source: [link](https://github.com/njoy/ACEtk) +- Dependencies: `pip install h5py numpy tqdm` +- You need a collection of ACE files (e.g., from NJOY or an ENDF/B distribution). ## Environment Variables | Variable | Description | |---------------|-------------------------------------------------------| -| `MCDC_ACELIB` | Path to the directory containing your ACE files. | -| `MCDC_LIB` | Path to the output directory for MC/DC HDF5 files. | +| `MCDC_ACELIB_NEUTRON` | Path to the directory containing your ACE files. | +| `MCDC_LIB_NEUTRON` | Path to the output directory for MC/DC HDF5 files. | ## Usage ```bash -export MCDC_ACELIB=/path/to/ace/files -export MCDC_LIB=/path/to/mcdc/library +export MCDC_ACELIB_NEUTRON=/path/to/ace/files +export MCDC_LIB_NEUTRON=/path/to/mcdc/library python generate.py # Convert only missing nuclides python generate.py --rewrite # Regenerate all files @@ -31,7 +29,7 @@ python generate.py --verbose # Print detailed per-nuclide info ## What it Does -For each ACE file in `$MCDC_ACELIB`, the generator: +For each ACE file in `$MCDC_ACELIB_NEUTRON`, the generator: 1. Reads the ACE header to identify the nuclide (Z, A, isomeric state) and temperature. 2. Extracts pointwise cross sections (elastic, capture, inelastic, fission) and the energy grid. diff --git a/tools/data_library_generator/generate.py b/tools/data_library_generator/neutron/generate.py similarity index 98% rename from tools/data_library_generator/generate.py rename to tools/data_library_generator/neutron/generate.py index 06a63dc76..94202f7fc 100644 --- a/tools/data_library_generator/generate.py +++ b/tools/data_library_generator/neutron/generate.py @@ -19,13 +19,13 @@ verbose = args.verbose # Directories -output_dir = os.getenv("MCDC_LIB") -ace_dir = os.getenv("MCDC_ACELIB") +output_dir = os.getenv("MCDC_LIB_NEUTRON") +ace_dir = os.getenv("MCDC_ACELIB_NEUTRON") if output_dir is None: - print_error("Environment variable $MCDC_LIB is not set") + print_error("Environment variable $MCDC_LIB_NEUTRON is not set") if ace_dir is None: - print_error("Environment variable $MCDC_ACELIB is not set") + print_error("Environment variable $MCDC_ACELIB_NEUTRON is not set") # Create output directory if needed os.makedirs(output_dir, exist_ok=True) diff --git a/tools/data_library_generator/util.py b/tools/data_library_generator/neutron/util.py similarity index 100% rename from tools/data_library_generator/util.py rename to tools/data_library_generator/neutron/util.py