diff --git a/mcdc/constant.py b/mcdc/constant.py index 79ade2223..2a8d20fa5 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -94,6 +94,7 @@ # Miscellanies EVENT_TIME_CENSUS = 1 << 5 EVENT_TIME_BOUNDARY = 1 << 6 +EVENT_CSDA_EDEP = 1 << 7 # Materials MATERIAL = 0 @@ -115,6 +116,10 @@ ELECTRON_REACTION_IONIZATION = 104 ELECTRON_REACTION_BREMSSTRAHLUNG = 105 ELECTRON_REACTION_EXCITATION = 106 +PROTON_REACTION_TOTAL = 200 +PROTON_REACTION_ELASTIC_SCATTERING = 201 +PROTON_REACTION_CAPTURE = 202 +PROTON_REACTION_INELASTIC_SCATTERING = 203 # Particle types PARTICLE_NEUTRON = 0 @@ -138,12 +143,12 @@ DISTRIBUTION_TABULATED_ENERGY_ANGLE = 8 DISTRIBUTION_N_BODY = 9 -# Anguler distribution type +# Angular distribution type ANGLE_ISOTROPIC = 0 ANGLE_DISTRIBUTED = 1 ANGLE_ENERGY_CORRELATED = 2 -# Referance frame +# Reference frame REFERENCE_FRAME_LAB = 0 REFERENCE_FRAME_COM = 1 @@ -193,8 +198,10 @@ LIGHT_SPEED = 2.99792458e10 # cm/s NEUTRON_MASS = 939.565413e6 # eV/c^2 ELECTRON_MASS = 510.99895069e3 # eV/c^2 +PROTON_MASS = 938.27208943e6 # eV/c^2 BOLTZMANN_K = 8.61733326e-5 # eV/K ELECTRON_CUTOFF_ENERGY = 100 # eV +PROTON_CUTOFF_ENERGY = 250000 # eV MU_CUTOFF = 0.999999 THERMAL_THRESHOLD_FACTOR = 400 diff --git a/mcdc/main.py b/mcdc/main.py index 71db344ae..3e621b467 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -168,6 +168,10 @@ def preparation(): if isinstance(material, Material): update_fissionable_from_nuclides(material) + if settings.proton_transport: + for nuclide in simulationPy.nuclides: + nuclide.set_proton_data() + if settings.electron_transport: for element in simulationPy.elements: element.set_electron_data() diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py index d1d03c34c..9a59d9c80 100644 --- a/mcdc/mcdc_get/__init__.py +++ b/mcdc/mcdc_get/__init__.py @@ -82,10 +82,18 @@ import mcdc.mcdc_get.neutron_inelastic_scattering_reaction as neutron_inelastic_scattering_reaction +import mcdc.mcdc_get.proton_capture_reaction as proton_capture_reaction + +import mcdc.mcdc_get.proton_elastic_scattering_reaction as proton_elastic_scattering_reaction + +import mcdc.mcdc_get.proton_inelastic_scattering_reaction as proton_inelastic_scattering_reaction + import mcdc.mcdc_get.collision_data as collision_data import mcdc.mcdc_get.particle_bank as particle_bank +import mcdc.mcdc_get.proton_reaction as proton_reaction + import mcdc.mcdc_get.settings as settings import mcdc.mcdc_get.global_weight_roulette as global_weight_roulette diff --git a/mcdc/mcdc_get/native_material.py b/mcdc/mcdc_get/native_material.py index 361294585..5aa48cbcb 100644 --- a/mcdc/mcdc_get/native_material.py +++ b/mcdc/mcdc_get/native_material.py @@ -117,3 +117,61 @@ def element_densities_chunk(start, length, native_material, data): start += native_material["element_densities_offset"] end = start + length return data[start:end] + + +@njit +def stopping_power(index, native_material, data): + offset = native_material["stopping_power_offset"] + return data[offset + index] + + +@njit +def stopping_power_all(native_material, data): + start = native_material["stopping_power_offset"] + size = native_material["stopping_power_length"] + end = start + size + return data[start:end] + + +@njit +def stopping_power_last(native_material, data): + start = native_material["stopping_power_offset"] + size = native_material["stopping_power_length"] + end = start + size + return data[end - 1] + + +@njit +def stopping_power_chunk(start, length, native_material, data): + start += native_material["stopping_power_offset"] + end = start + length + return data[start:end] + + +@njit +def stopping_power_energy_grid(index, native_material, data): + offset = native_material["stopping_power_energy_grid_offset"] + return data[offset + index] + + +@njit +def stopping_power_energy_grid_all(native_material, data): + start = native_material["stopping_power_energy_grid_offset"] + size = native_material["stopping_power_energy_grid_length"] + end = start + size + return data[start:end] + + +@njit +def stopping_power_energy_grid_last(native_material, data): + start = native_material["stopping_power_energy_grid_offset"] + size = native_material["stopping_power_energy_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def stopping_power_energy_grid_chunk(start, length, native_material, data): + start += native_material["stopping_power_energy_grid_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/nuclide.py b/mcdc/mcdc_get/nuclide.py index af9593f10..e30d775dd 100644 --- a/mcdc/mcdc_get/nuclide.py +++ b/mcdc/mcdc_get/nuclide.py @@ -177,6 +177,151 @@ def neutron_fission_xs_chunk(start, length, nuclide, data): return data[start:end] +@njit +def proton_xs_energy_grid(index, nuclide, data): + offset = nuclide["proton_xs_energy_grid_offset"] + return data[offset + index] + + +@njit +def proton_xs_energy_grid_all(nuclide, data): + start = nuclide["proton_xs_energy_grid_offset"] + size = nuclide["proton_xs_energy_grid_length"] + end = start + size + return data[start:end] + + +@njit +def proton_xs_energy_grid_last(nuclide, data): + start = nuclide["proton_xs_energy_grid_offset"] + size = nuclide["proton_xs_energy_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def proton_xs_energy_grid_chunk(start, length, nuclide, data): + start += nuclide["proton_xs_energy_grid_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_total_xs(index, nuclide, data): + offset = nuclide["proton_total_xs_offset"] + return data[offset + index] + + +@njit +def proton_total_xs_all(nuclide, data): + start = nuclide["proton_total_xs_offset"] + size = nuclide["proton_total_xs_length"] + end = start + size + return data[start:end] + + +@njit +def proton_total_xs_last(nuclide, data): + start = nuclide["proton_total_xs_offset"] + size = nuclide["proton_total_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def proton_total_xs_chunk(start, length, nuclide, data): + start += nuclide["proton_total_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_elastic_xs(index, nuclide, data): + offset = nuclide["proton_elastic_xs_offset"] + return data[offset + index] + + +@njit +def proton_elastic_xs_all(nuclide, data): + start = nuclide["proton_elastic_xs_offset"] + size = nuclide["proton_elastic_xs_length"] + end = start + size + return data[start:end] + + +@njit +def proton_elastic_xs_last(nuclide, data): + start = nuclide["proton_elastic_xs_offset"] + size = nuclide["proton_elastic_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def proton_elastic_xs_chunk(start, length, nuclide, data): + start += nuclide["proton_elastic_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_capture_xs(index, nuclide, data): + offset = nuclide["proton_capture_xs_offset"] + return data[offset + index] + + +@njit +def proton_capture_xs_all(nuclide, data): + start = nuclide["proton_capture_xs_offset"] + size = nuclide["proton_capture_xs_length"] + end = start + size + return data[start:end] + + +@njit +def proton_capture_xs_last(nuclide, data): + start = nuclide["proton_capture_xs_offset"] + size = nuclide["proton_capture_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def proton_capture_xs_chunk(start, length, nuclide, data): + start += nuclide["proton_capture_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_inelastic_xs(index, nuclide, data): + offset = nuclide["proton_inelastic_xs_offset"] + return data[offset + index] + + +@njit +def proton_inelastic_xs_all(nuclide, data): + start = nuclide["proton_inelastic_xs_offset"] + size = nuclide["proton_inelastic_xs_length"] + end = start + size + return data[start:end] + + +@njit +def proton_inelastic_xs_last(nuclide, data): + start = nuclide["proton_inelastic_xs_offset"] + size = nuclide["proton_inelastic_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def proton_inelastic_xs_chunk(start, length, nuclide, data): + start += nuclide["proton_inelastic_xs_offset"] + end = start + length + return data[start:end] + + @njit def neutron_elastic_scattering_reaction_IDs(index, nuclide, data): offset = nuclide["neutron_elastic_scattering_reaction_IDs_offset"] @@ -293,6 +438,93 @@ def neutron_fission_reaction_IDs_chunk(start, length, nuclide, data): return data[start:end] +@njit +def proton_elastic_scattering_reaction_IDs(index, nuclide, data): + offset = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def proton_elastic_scattering_reaction_IDs_all(nuclide, data): + start = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_elastic_scattering_reaction"] + end = start + size + return data[start:end] + + +@njit +def proton_elastic_scattering_reaction_IDs_last(nuclide, data): + start = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_elastic_scattering_reaction"] + end = start + size + return data[end - 1] + + +@njit +def proton_elastic_scattering_reaction_IDs_chunk(start, length, nuclide, data): + start += nuclide["proton_elastic_scattering_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_capture_reaction_IDs(index, nuclide, data): + offset = nuclide["proton_capture_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def proton_capture_reaction_IDs_all(nuclide, data): + start = nuclide["proton_capture_reaction_IDs_offset"] + size = nuclide["N_proton_capture_reaction"] + end = start + size + return data[start:end] + + +@njit +def proton_capture_reaction_IDs_last(nuclide, data): + start = nuclide["proton_capture_reaction_IDs_offset"] + size = nuclide["N_proton_capture_reaction"] + end = start + size + return data[end - 1] + + +@njit +def proton_capture_reaction_IDs_chunk(start, length, nuclide, data): + start += nuclide["proton_capture_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def proton_inelastic_scattering_reaction_IDs(index, nuclide, data): + offset = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def proton_inelastic_scattering_reaction_IDs_all(nuclide, data): + start = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_inelastic_scattering_reaction"] + end = start + size + return data[start:end] + + +@njit +def proton_inelastic_scattering_reaction_IDs_last(nuclide, data): + start = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_inelastic_scattering_reaction"] + end = start + size + return data[end - 1] + + +@njit +def proton_inelastic_scattering_reaction_IDs_chunk(start, length, nuclide, data): + start += nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + end = start + length + return data[start:end] + + @njit def neutron_fission_delayed_fractions(index, nuclide, data): offset = nuclide["neutron_fission_delayed_fractions_offset"] @@ -378,3 +610,61 @@ def neutron_fission_delayed_spectrum_IDs_chunk(start, length, nuclide, data): start += nuclide["neutron_fission_delayed_spectrum_IDs_offset"] end = start + length return data[start:end] + + +@njit +def stopping_power(index, nuclide, data): + offset = nuclide["stopping_power_offset"] + return data[offset + index] + + +@njit +def stopping_power_all(nuclide, data): + start = nuclide["stopping_power_offset"] + size = nuclide["stopping_power_length"] + end = start + size + return data[start:end] + + +@njit +def stopping_power_last(nuclide, data): + start = nuclide["stopping_power_offset"] + size = nuclide["stopping_power_length"] + end = start + size + return data[end - 1] + + +@njit +def stopping_power_chunk(start, length, nuclide, data): + start += nuclide["stopping_power_offset"] + end = start + length + return data[start:end] + + +@njit +def stopping_power_energy_grid(index, nuclide, data): + offset = nuclide["stopping_power_energy_grid_offset"] + return data[offset + index] + + +@njit +def stopping_power_energy_grid_all(nuclide, data): + start = nuclide["stopping_power_energy_grid_offset"] + size = nuclide["stopping_power_energy_grid_length"] + end = start + size + return data[start:end] + + +@njit +def stopping_power_energy_grid_last(nuclide, data): + start = nuclide["stopping_power_energy_grid_offset"] + size = nuclide["stopping_power_energy_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def stopping_power_energy_grid_chunk(start, length, nuclide, data): + start += nuclide["stopping_power_energy_grid_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/proton_capture_reaction.py b/mcdc/mcdc_get/proton_capture_reaction.py new file mode 100644 index 000000000..fdbf8e750 --- /dev/null +++ b/mcdc/mcdc_get/proton_capture_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/proton_elastic_scattering_reaction.py b/mcdc/mcdc_get/proton_elastic_scattering_reaction.py new file mode 100644 index 000000000..fdbf8e750 --- /dev/null +++ b/mcdc/mcdc_get/proton_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/proton_inelastic_scattering_reaction.py b/mcdc/mcdc_get/proton_inelastic_scattering_reaction.py new file mode 100644 index 000000000..abb8e860a --- /dev/null +++ b/mcdc/mcdc_get/proton_inelastic_scattering_reaction.py @@ -0,0 +1,84 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def spectrum_probability_grid(index, proton_inelastic_scattering_reaction, data): + offset = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + return data[offset + index] + + +@njit +def spectrum_probability_grid_all(proton_inelastic_scattering_reaction, data): + start = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + size = proton_inelastic_scattering_reaction["spectrum_probability_grid_length"] + end = start + size + return data[start:end] + + +@njit +def spectrum_probability_grid_last(proton_inelastic_scattering_reaction, data): + start = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + size = proton_inelastic_scattering_reaction["spectrum_probability_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def spectrum_probability_grid_chunk(start, length, proton_inelastic_scattering_reaction, data): + start += proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + end = start + length + return data[start:end] + + +@njit +def spectrum_probability_vector(index_1, proton_inelastic_scattering_reaction, data): + offset = proton_inelastic_scattering_reaction["spectrum_probability_offset"] + stride = proton_inelastic_scattering_reaction["N_spectrum"] + start = offset + index_1 * stride + end = start + stride + return data[start:end] + + +@njit +def spectrum_probability(index_1, index_2, proton_inelastic_scattering_reaction, data): + offset = proton_inelastic_scattering_reaction["spectrum_probability_offset"] + stride = proton_inelastic_scattering_reaction["N_spectrum"] + return data[offset + index_1 * stride + index_2] + + +@njit +def spectrum_probability_chunk(start, length, proton_inelastic_scattering_reaction, data): + start += proton_inelastic_scattering_reaction["spectrum_probability_offset"] + end = start + length + return data[start:end] + + +@njit +def energy_spectrum_IDs(index, proton_inelastic_scattering_reaction, data): + offset = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + return data[offset + index] + + +@njit +def energy_spectrum_IDs_all(proton_inelastic_scattering_reaction, data): + start = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + size = proton_inelastic_scattering_reaction["N_energy_spectrum"] + end = start + size + return data[start:end] + + +@njit +def energy_spectrum_IDs_last(proton_inelastic_scattering_reaction, data): + start = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + size = proton_inelastic_scattering_reaction["N_energy_spectrum"] + end = start + size + return data[end - 1] + + +@njit +def energy_spectrum_IDs_chunk(start, length, proton_inelastic_scattering_reaction, data): + start += proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/proton_reaction.py b/mcdc/mcdc_get/proton_reaction.py new file mode 100644 index 000000000..94fd03787 --- /dev/null +++ b/mcdc/mcdc_get/proton_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs(index, proton_reaction, data): + offset = proton_reaction["xs_offset"] + return data[offset + index] + + +@njit +def xs_all(proton_reaction, data): + start = proton_reaction["xs_offset"] + size = proton_reaction["xs_length"] + end = start + size + return data[start:end] + + +@njit +def xs_last(proton_reaction, data): + start = proton_reaction["xs_offset"] + size = proton_reaction["xs_length"] + end = start + size + return data[end - 1] + + +@njit +def xs_chunk(start, length, proton_reaction, data): + start += proton_reaction["xs_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py index 8f771c058..8395de1e6 100644 --- a/mcdc/mcdc_set/__init__.py +++ b/mcdc/mcdc_set/__init__.py @@ -82,10 +82,18 @@ import mcdc.mcdc_set.neutron_inelastic_scattering_reaction as neutron_inelastic_scattering_reaction +import mcdc.mcdc_set.proton_capture_reaction as proton_capture_reaction + +import mcdc.mcdc_set.proton_elastic_scattering_reaction as proton_elastic_scattering_reaction + +import mcdc.mcdc_set.proton_inelastic_scattering_reaction as proton_inelastic_scattering_reaction + import mcdc.mcdc_set.collision_data as collision_data import mcdc.mcdc_set.particle_bank as particle_bank +import mcdc.mcdc_set.proton_reaction as proton_reaction + import mcdc.mcdc_set.settings as settings import mcdc.mcdc_set.global_weight_roulette as global_weight_roulette diff --git a/mcdc/mcdc_set/native_material.py b/mcdc/mcdc_set/native_material.py index 303ffb474..ea7845d84 100644 --- a/mcdc/mcdc_set/native_material.py +++ b/mcdc/mcdc_set/native_material.py @@ -117,3 +117,61 @@ def element_densities_chunk(start, length, native_material, data, value): start += native_material["element_densities_offset"] end = start + length data[start:end] = value + + +@njit +def stopping_power(index, native_material, data, value): + offset = native_material["stopping_power_offset"] + data[offset + index] = value + + +@njit +def stopping_power_all(native_material, data, value): + start = native_material["stopping_power_offset"] + size = native_material["stopping_power_length"] + end = start + size + data[start:end] = value + + +@njit +def stopping_power_last(native_material, data, value): + start = native_material["stopping_power_offset"] + size = native_material["stopping_power_length"] + end = start + size + data[end - 1] = value + + +@njit +def stopping_power_chunk(start, length, native_material, data, value): + start += native_material["stopping_power_offset"] + end = start + length + data[start:end] = value + + +@njit +def stopping_power_energy_grid(index, native_material, data, value): + offset = native_material["stopping_power_energy_grid_offset"] + data[offset + index] = value + + +@njit +def stopping_power_energy_grid_all(native_material, data, value): + start = native_material["stopping_power_energy_grid_offset"] + size = native_material["stopping_power_energy_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def stopping_power_energy_grid_last(native_material, data, value): + start = native_material["stopping_power_energy_grid_offset"] + size = native_material["stopping_power_energy_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def stopping_power_energy_grid_chunk(start, length, native_material, data, value): + start += native_material["stopping_power_energy_grid_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/nuclide.py b/mcdc/mcdc_set/nuclide.py index 257d62580..c9d60e0d2 100644 --- a/mcdc/mcdc_set/nuclide.py +++ b/mcdc/mcdc_set/nuclide.py @@ -177,6 +177,151 @@ def neutron_fission_xs_chunk(start, length, nuclide, data, value): data[start:end] = value +@njit +def proton_xs_energy_grid(index, nuclide, data, value): + offset = nuclide["proton_xs_energy_grid_offset"] + data[offset + index] = value + + +@njit +def proton_xs_energy_grid_all(nuclide, data, value): + start = nuclide["proton_xs_energy_grid_offset"] + size = nuclide["proton_xs_energy_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def proton_xs_energy_grid_last(nuclide, data, value): + start = nuclide["proton_xs_energy_grid_offset"] + size = nuclide["proton_xs_energy_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def proton_xs_energy_grid_chunk(start, length, nuclide, data, value): + start += nuclide["proton_xs_energy_grid_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_total_xs(index, nuclide, data, value): + offset = nuclide["proton_total_xs_offset"] + data[offset + index] = value + + +@njit +def proton_total_xs_all(nuclide, data, value): + start = nuclide["proton_total_xs_offset"] + size = nuclide["proton_total_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def proton_total_xs_last(nuclide, data, value): + start = nuclide["proton_total_xs_offset"] + size = nuclide["proton_total_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def proton_total_xs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_total_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_elastic_xs(index, nuclide, data, value): + offset = nuclide["proton_elastic_xs_offset"] + data[offset + index] = value + + +@njit +def proton_elastic_xs_all(nuclide, data, value): + start = nuclide["proton_elastic_xs_offset"] + size = nuclide["proton_elastic_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def proton_elastic_xs_last(nuclide, data, value): + start = nuclide["proton_elastic_xs_offset"] + size = nuclide["proton_elastic_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def proton_elastic_xs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_elastic_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_capture_xs(index, nuclide, data, value): + offset = nuclide["proton_capture_xs_offset"] + data[offset + index] = value + + +@njit +def proton_capture_xs_all(nuclide, data, value): + start = nuclide["proton_capture_xs_offset"] + size = nuclide["proton_capture_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def proton_capture_xs_last(nuclide, data, value): + start = nuclide["proton_capture_xs_offset"] + size = nuclide["proton_capture_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def proton_capture_xs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_capture_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_inelastic_xs(index, nuclide, data, value): + offset = nuclide["proton_inelastic_xs_offset"] + data[offset + index] = value + + +@njit +def proton_inelastic_xs_all(nuclide, data, value): + start = nuclide["proton_inelastic_xs_offset"] + size = nuclide["proton_inelastic_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def proton_inelastic_xs_last(nuclide, data, value): + start = nuclide["proton_inelastic_xs_offset"] + size = nuclide["proton_inelastic_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def proton_inelastic_xs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_inelastic_xs_offset"] + end = start + length + data[start:end] = value + + @njit def neutron_elastic_scattering_reaction_IDs(index, nuclide, data, value): offset = nuclide["neutron_elastic_scattering_reaction_IDs_offset"] @@ -293,6 +438,93 @@ def neutron_fission_reaction_IDs_chunk(start, length, nuclide, data, value): data[start:end] = value +@njit +def proton_elastic_scattering_reaction_IDs(index, nuclide, data, value): + offset = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def proton_elastic_scattering_reaction_IDs_all(nuclide, data, value): + start = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_elastic_scattering_reaction"] + end = start + size + data[start:end] = value + + +@njit +def proton_elastic_scattering_reaction_IDs_last(nuclide, data, value): + start = nuclide["proton_elastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_elastic_scattering_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def proton_elastic_scattering_reaction_IDs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_elastic_scattering_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_capture_reaction_IDs(index, nuclide, data, value): + offset = nuclide["proton_capture_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def proton_capture_reaction_IDs_all(nuclide, data, value): + start = nuclide["proton_capture_reaction_IDs_offset"] + size = nuclide["N_proton_capture_reaction"] + end = start + size + data[start:end] = value + + +@njit +def proton_capture_reaction_IDs_last(nuclide, data, value): + start = nuclide["proton_capture_reaction_IDs_offset"] + size = nuclide["N_proton_capture_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def proton_capture_reaction_IDs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_capture_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def proton_inelastic_scattering_reaction_IDs(index, nuclide, data, value): + offset = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def proton_inelastic_scattering_reaction_IDs_all(nuclide, data, value): + start = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_inelastic_scattering_reaction"] + end = start + size + data[start:end] = value + + +@njit +def proton_inelastic_scattering_reaction_IDs_last(nuclide, data, value): + start = nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + size = nuclide["N_proton_inelastic_scattering_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def proton_inelastic_scattering_reaction_IDs_chunk(start, length, nuclide, data, value): + start += nuclide["proton_inelastic_scattering_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + @njit def neutron_fission_delayed_fractions(index, nuclide, data, value): offset = nuclide["neutron_fission_delayed_fractions_offset"] @@ -378,3 +610,61 @@ def neutron_fission_delayed_spectrum_IDs_chunk(start, length, nuclide, data, val start += nuclide["neutron_fission_delayed_spectrum_IDs_offset"] end = start + length data[start:end] = value + + +@njit +def stopping_power(index, nuclide, data, value): + offset = nuclide["stopping_power_offset"] + data[offset + index] = value + + +@njit +def stopping_power_all(nuclide, data, value): + start = nuclide["stopping_power_offset"] + size = nuclide["stopping_power_length"] + end = start + size + data[start:end] = value + + +@njit +def stopping_power_last(nuclide, data, value): + start = nuclide["stopping_power_offset"] + size = nuclide["stopping_power_length"] + end = start + size + data[end - 1] = value + + +@njit +def stopping_power_chunk(start, length, nuclide, data, value): + start += nuclide["stopping_power_offset"] + end = start + length + data[start:end] = value + + +@njit +def stopping_power_energy_grid(index, nuclide, data, value): + offset = nuclide["stopping_power_energy_grid_offset"] + data[offset + index] = value + + +@njit +def stopping_power_energy_grid_all(nuclide, data, value): + start = nuclide["stopping_power_energy_grid_offset"] + size = nuclide["stopping_power_energy_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def stopping_power_energy_grid_last(nuclide, data, value): + start = nuclide["stopping_power_energy_grid_offset"] + size = nuclide["stopping_power_energy_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def stopping_power_energy_grid_chunk(start, length, nuclide, data, value): + start += nuclide["stopping_power_energy_grid_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/proton_capture_reaction.py b/mcdc/mcdc_set/proton_capture_reaction.py new file mode 100644 index 000000000..fdbf8e750 --- /dev/null +++ b/mcdc/mcdc_set/proton_capture_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/proton_elastic_scattering_reaction.py b/mcdc/mcdc_set/proton_elastic_scattering_reaction.py new file mode 100644 index 000000000..fdbf8e750 --- /dev/null +++ b/mcdc/mcdc_set/proton_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/proton_inelastic_scattering_reaction.py b/mcdc/mcdc_set/proton_inelastic_scattering_reaction.py new file mode 100644 index 000000000..485e20529 --- /dev/null +++ b/mcdc/mcdc_set/proton_inelastic_scattering_reaction.py @@ -0,0 +1,84 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def spectrum_probability_grid(index, proton_inelastic_scattering_reaction, data, value): + offset = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + data[offset + index] = value + + +@njit +def spectrum_probability_grid_all(proton_inelastic_scattering_reaction, data, value): + start = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + size = proton_inelastic_scattering_reaction["spectrum_probability_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def spectrum_probability_grid_last(proton_inelastic_scattering_reaction, data, value): + start = proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + size = proton_inelastic_scattering_reaction["spectrum_probability_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def spectrum_probability_grid_chunk(start, length, proton_inelastic_scattering_reaction, data, value): + start += proton_inelastic_scattering_reaction["spectrum_probability_grid_offset"] + end = start + length + data[start:end] = value + + +@njit +def spectrum_probability_vector(index_1, proton_inelastic_scattering_reaction, data, value): + offset = proton_inelastic_scattering_reaction["spectrum_probability_offset"] + stride = proton_inelastic_scattering_reaction["N_spectrum"] + start = offset + index_1 * stride + end = start + stride + data[start:end] = value + + +@njit +def spectrum_probability(index_1, index_2, proton_inelastic_scattering_reaction, data, value): + offset = proton_inelastic_scattering_reaction["spectrum_probability_offset"] + stride = proton_inelastic_scattering_reaction["N_spectrum"] + data[offset + index_1 * stride + index_2] = value + + +@njit +def spectrum_probability_chunk(start, length, proton_inelastic_scattering_reaction, data, value): + start += proton_inelastic_scattering_reaction["spectrum_probability_offset"] + end = start + length + data[start:end] = value + + +@njit +def energy_spectrum_IDs(index, proton_inelastic_scattering_reaction, data, value): + offset = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + data[offset + index] = value + + +@njit +def energy_spectrum_IDs_all(proton_inelastic_scattering_reaction, data, value): + start = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + size = proton_inelastic_scattering_reaction["N_energy_spectrum"] + end = start + size + data[start:end] = value + + +@njit +def energy_spectrum_IDs_last(proton_inelastic_scattering_reaction, data, value): + start = proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + size = proton_inelastic_scattering_reaction["N_energy_spectrum"] + end = start + size + data[end - 1] = value + + +@njit +def energy_spectrum_IDs_chunk(start, length, proton_inelastic_scattering_reaction, data, value): + start += proton_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/proton_reaction.py b/mcdc/mcdc_set/proton_reaction.py new file mode 100644 index 000000000..5221e64d5 --- /dev/null +++ b/mcdc/mcdc_set/proton_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs(index, proton_reaction, data, value): + offset = proton_reaction["xs_offset"] + data[offset + index] = value + + +@njit +def xs_all(proton_reaction, data, value): + start = proton_reaction["xs_offset"] + size = proton_reaction["xs_length"] + end = start + size + data[start:end] = value + + +@njit +def xs_last(proton_reaction, data, value): + start = proton_reaction["xs_offset"] + size = proton_reaction["xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def xs_chunk(start, length, proton_reaction, data, value): + start += proton_reaction["xs_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py index 0948870c6..52fe42146 100644 --- a/mcdc/numba_types.py +++ b/mcdc/numba_types.py @@ -342,6 +342,13 @@ ('nuclide_densities_length', int64), ('element_densities_offset', int64), ('element_densities_length', int64), + ('stopping_power_provided', bool), + ('stopping_power_offset', int64), + ('stopping_power_length', int64), + ('stopping_power_energy_grid_offset', int64), + ('stopping_power_energy_grid_length', int64), + ('radiation_length', float64), + ('radiation_length_provided', bool), ('ID', int64), ('parent_ID', int64), ]) @@ -389,6 +396,7 @@ ('atomic_weight_ratio', float64), ('fissionable', bool), ('excitation_level', int64), + ('radiation_length', float64), ('neutron_xs_energy_grid_offset', int64), ('neutron_xs_energy_grid_length', int64), ('neutron_total_xs_offset', int64), @@ -401,6 +409,16 @@ ('neutron_inelastic_xs_length', int64), ('neutron_fission_xs_offset', int64), ('neutron_fission_xs_length', int64), + ('proton_xs_energy_grid_offset', int64), + ('proton_xs_energy_grid_length', int64), + ('proton_total_xs_offset', int64), + ('proton_total_xs_length', int64), + ('proton_elastic_xs_offset', int64), + ('proton_elastic_xs_length', int64), + ('proton_capture_xs_offset', int64), + ('proton_capture_xs_length', int64), + ('proton_inelastic_xs_offset', int64), + ('proton_inelastic_xs_length', int64), ('N_neutron_elastic_scattering_reaction', int64), ('neutron_elastic_scattering_reaction_IDs_offset', int64), ('N_neutron_capture_reaction', int64), @@ -409,6 +427,12 @@ ('neutron_inelastic_scattering_reaction_IDs_offset', int64), ('N_neutron_fission_reaction', int64), ('neutron_fission_reaction_IDs_offset', int64), + ('N_proton_elastic_scattering_reaction', int64), + ('proton_elastic_scattering_reaction_IDs_offset', int64), + ('N_proton_capture_reaction', int64), + ('proton_capture_reaction_IDs_offset', int64), + ('N_proton_inelastic_scattering_reaction', int64), + ('proton_inelastic_scattering_reaction_IDs_offset', int64), ('neutron_fission_prompt_multiplicity_ID', int64), ('neutron_fission_delayed_multiplicity_ID', int64), ('N_neutron_fission_delayed_precursor', int64), @@ -418,6 +442,10 @@ ('neutron_fission_delayed_decay_rates_length', int64), ('N_neutron_fission_delayed_spectrum', int64), ('neutron_fission_delayed_spectrum_IDs_offset', int64), + ('stopping_power_offset', int64), + ('stopping_power_length', int64), + ('stopping_power_energy_grid_offset', int64), + ('stopping_power_energy_grid_length', int64), ('ID', int64), ]) @@ -504,6 +532,33 @@ ('parent_ID', int64), ]) +proton_capture_reaction = into_dtype([ + ('ID', int64), + ('parent_ID', int64), +]) + +proton_elastic_scattering_reaction = into_dtype([ + ('mu_table_ID', int64), + ('ID', int64), + ('parent_ID', int64), +]) + +proton_inelastic_scattering_reaction = into_dtype([ + ('multiplicity', int64), + ('angle_type', int64), + ('mu_ID', int64), + ('N_spectrum_probability_bin', int64), + ('N_spectrum', int64), + ('spectrum_probability_grid_offset', int64), + ('spectrum_probability_grid_length', int64), + ('spectrum_probability_offset', int64), + ('spectrum_probability_length', int64), + ('N_energy_spectrum', int64), + ('energy_spectrum_IDs_offset', int64), + ('ID', int64), + ('parent_ID', int64), +]) + collision_data = into_dtype([ ('energy_deposition', float64), ]) @@ -513,6 +568,18 @@ ('tag', 'U32'), ]) +proton_reaction = into_dtype([ + ('MT', int64), + ('xs_offset', int64), + ('xs_length', int64), + ('xs_offset_', int64), + ('reference_frame', int64), + ('q_value', float64), + ('ID', int64), + ('child_type', int64), + ('child_ID', int64), +]) + settings = into_dtype([ ('N_particle', int64), ('N_batch', int64), @@ -528,6 +595,8 @@ ('time_boundary', float64), ('output_name', 'U32'), ('use_progress_bar', bool), + ('csda', bool), + ('csda_max_fractional_e_loss', float64), ('N_census', int64), ('census_time_offset', int64), ('census_time_length', int64), @@ -810,6 +879,14 @@ def set_simulation(N: dict): ('N_neutron_inelastic_scattering_reaction', int64), ('sources', source, (N['source'])), ('N_source', int64), + ('proton_capture_reactions', proton_capture_reaction, (N['proton_capture_reaction'])), + ('N_proton_capture_reaction', int64), + ('proton_elastic_scattering_reactions', proton_elastic_scattering_reaction, (N['proton_elastic_scattering_reaction'])), + ('N_proton_elastic_scattering_reaction', int64), + ('proton_inelastic_scattering_reactions', proton_inelastic_scattering_reaction, (N['proton_inelastic_scattering_reaction'])), + ('N_proton_inelastic_scattering_reaction', int64), + ('proton_reactions', proton_reaction, (N['proton_reaction'])), + ('N_proton_reaction', int64), ('cells', cell, (N['cell'])), ('N_cell', int64), ('lattices', lattice, (N['lattice'])), diff --git a/mcdc/object_/base.py b/mcdc/object_/base.py index 11d87e8ad..bc52e3a05 100644 --- a/mcdc/object_/base.py +++ b/mcdc/object_/base.py @@ -70,6 +70,7 @@ def register_object(object_): from mcdc.object_.mesh import MeshBase from mcdc.object_.nuclide import Nuclide from mcdc.object_.neutron_reaction import NeutronReactionBase + from mcdc.object_.proton_reaction import ProtonReactionBase from mcdc.object_.source import Source from mcdc.object_.surface import Surface from mcdc.object_.tally import Tally @@ -95,6 +96,8 @@ def register_object(object_): object_list = simulation.nuclides elif isinstance(object_, NeutronReactionBase): object_list = simulation.neutron_reactions + elif isinstance(object_, ProtonReactionBase): + object_list = simulation.proton_reactions elif isinstance(object_, Region): object_list = simulation.regions elif isinstance(object_, Source): diff --git a/mcdc/object_/material.py b/mcdc/object_/material.py index c2b2fbfbe..e1510b5ce 100644 --- a/mcdc/object_/material.py +++ b/mcdc/object_/material.py @@ -1,5 +1,7 @@ import numpy as np import os +import h5py +import re from numpy import float64 from numpy.typing import NDArray @@ -69,9 +71,9 @@ class Material(MaterialBase): name : str, optional User label. nuclide_composition : dict - Dictionary mapping nuclide names (str) to atom densities (float). + Dictionary mapping nuclide names (str) to atom densities in units of atoms/barn-cm (float). element_composition : dict - Dictionary mapping element names (str) to atom densities (float). + Dictionary mapping element names (str) to atom densities in units of atoms/barn-cm (float). temperature : float, optional Temperature in Kelvin (default 293.6 K). @@ -101,6 +103,13 @@ class Material(MaterialBase): elements: list[Element] nuclide_densities: NDArray[float64] element_densities: NDArray[float64] + # + stopping_power_provided: bool = False + stopping_power: NDArray[float64] + stopping_power_energy_grid: NDArray[float64] + # + radiation_length: float = 0.0 + radiation_length_provided: bool = False def __init__( self, @@ -129,6 +138,14 @@ def __init__( self.elements = [] self.element_densities = np.zeros(len(element_composition)) + # Stopping power + self.stopping_power = np.array([]) + self.stopping_power_energy_grid = np.array([]) + + # Radiation length calculation prep + total_mass = 0.0 + X0_weighted_mass = 0.0 + # Check if library directory is set lib_dir = os.getenv("MCDC_LIB") if lib_dir is None: @@ -204,6 +221,16 @@ def __init__( if nuclide.fissionable: self.fissionable = True + # Calculate the material's radiation length (for proton transport purposes) + nuclide_mass = nuclide.mass_number + nuclide_X0 = nuclide.radiation_length + + total_mass += nuclide_mass * nuclide_density + X0_weighted_mass += nuclide_mass * nuclide_density / nuclide_X0 + + # Set the material radiation length + self.radiation_length = total_mass / X0_weighted_mass + def __repr__(self): text = super().__repr__() text += f" - Temperature: {self.temperature} K\n" @@ -220,10 +247,35 @@ def __repr__(self): f" - {element.name:<5} | {self.element_composition[element]}\n" ) return text + + def add_stopping_power( + self, + stopping_power_filename: str = "", + ): + + self.stopping_power_provided = True + + dir_name = os.getenv("MCDC_LIB") + file_name = stopping_power_filename + file = h5py.File(f"{dir_name}/{file_name}.h5", "r") + + self.stopping_power = file["stopping_power"]["total_stopping_power"][()] + self.stopping_power_energy_grid = file["stopping_power"]["energy"][()] + if file["radiation_length"]["radiation_length"][()]: + self.radiation_length = file["radiation_length"]["radiation_length"][()] + file.close() + + def custom_radiation_length( + self, + radiation_length: float, + ): + + self.radiation_length_provided = True + self.radiation_length = radiation_length # Currently supported temperatures -TEMPERATURES = [0.1, 233.15, 273.15, 293.6, 600.0, 900.0, 1200.0, 2500.0] +TEMPERATURES = [0.0, 0.1, 233.15, 273.15, 293.6, 600.0, 900.0, 1200.0, 2500.0] # ====================================================================================== @@ -447,6 +499,7 @@ def __repr__(self): return text + def set_nuclides_from_elements(material): material.nuclides = [] material.nuclide_composition = {} diff --git a/mcdc/object_/nuclide.py b/mcdc/object_/nuclide.py index 80ca3e669..0c82ceef2 100644 --- a/mcdc/object_/nuclide.py +++ b/mcdc/object_/nuclide.py @@ -18,6 +18,12 @@ NeutronReactionInelasticScattering, set_energy_distribution, ) +from mcdc.object_.proton_reaction import ( + ProtonReactionElasticScattering, + ProtonReactionInelasticScattering, + ProtonReactionCapture, + set_energy_distribution, +) from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error @@ -37,6 +43,7 @@ class Nuclide(ObjectNonSingleton): atomic_weight_ratio: float fissionable: bool excitation_level: int + radiation_length: float # neutron_xs_energy_grid: NDArray[float64] neutron_total_xs: NDArray[float64] @@ -45,17 +52,30 @@ class Nuclide(ObjectNonSingleton): neutron_inelastic_xs: NDArray[float64] neutron_fission_xs: NDArray[float64] # + proton_xs_energy_grid: NDArray[float64] + proton_total_xs: NDArray[float64] + proton_elastic_xs: NDArray[float64] + proton_capture_xs: NDArray[float64] + proton_inelastic_xs: NDArray[float64] + # neutron_elastic_scattering_reactions: list[NeutronReactionElasticScattering] neutron_capture_reactions: list[NeutronReactionCapture] neutron_inelastic_scattering_reactions: list[NeutronReactionInelasticScattering] neutron_fission_reactions: list[NeutronReactionFission] # + proton_elastic_scattering_reactions: list[ProtonReactionElasticScattering] + proton_capture_reactions: list[ProtonReactionCapture] + proton_inelastic_scattering_reactions: list[ProtonReactionInelasticScattering] + # neutron_fission_prompt_multiplicity: DataBase neutron_fission_delayed_multiplicity: DataBase N_neutron_fission_delayed_precursor: int neutron_fission_delayed_fractions: NDArray[float64] neutron_fission_delayed_decay_rates: NDArray[float64] neutron_fission_delayed_spectra: list[DistributionBase] + # + stopping_power: NDArray[float64] + stopping_power_energy_grid: NDArray[float64] def __init__(self, nuclide_name, temperature): super().__init__() @@ -72,8 +92,42 @@ def __init__(self, nuclide_name, temperature): self.atomic_weight_ratio = file["atomic_weight_ratio"][()] self.fissionable = bool(file["fissionable"][()]) self.excitation_level = int(file["excitation_level"][()]) + self.radiation_length = float(file["radiation_length"][()]) file.close() + # Initialize all attributes to defaults + # Neutron XS + self.neutron_xs_energy_grid = np.zeros(0) + self.neutron_total_xs = np.zeros(0) + self.neutron_elastic_xs = np.zeros(0) + self.neutron_capture_xs = np.zeros(0) + self.neutron_inelastic_xs = np.zeros(0) + self.neutron_fission_xs = np.zeros(0) + # Proton XS + self.proton_xs_energy_grid = np.zeros(0) + self.proton_total_xs = np.zeros(0) + self.proton_elastic_xs = np.zeros(0) + self.proton_inelastic_xs = np.zeros(0) + self.proton_capture_xs = np.zeros(0) + # Reactions + self.neutron_elastic_scattering_reactions = [] + self.neutron_capture_reactions = [] + self.neutron_inelastic_scattering_reactions = [] + self.neutron_fission_reactions = [] + self.proton_elastic_scattering_reactions = [] + self.proton_inelastic_scattering_reactions = [] + self.proton_capture_reactions = [] + # Fission + self.neutron_fission_prompt_multiplicity = DataPolynomial(np.array([0.0])) + self.neutron_fission_delayed_multiplicity = DataPolynomial(np.array([0.0])) + self.N_neutron_fission_delayed_precursor = 0 + self.neutron_fission_delayed_fractions = np.zeros(0) + self.neutron_fission_delayed_decay_rates = np.zeros(0) + self.neutron_fission_delayed_spectra = [] + # Stopping Power + self.stopping_power = np.zeros(0) + self.stopping_power_energy_grid = np.zeros(0) + def set_neutron_data(self): nuclide_name = self.name temperature = self.temperature @@ -213,6 +267,110 @@ def set_neutron_data(self): file.close() + def set_proton_data(self): + nuclide_name = self.name + temperature = self.temperature + + # Load data library + dir_name = os.getenv("MCDC_LIB") + file_name = f"{nuclide_name}-{temperature}K.h5" + file = h5py.File(f"{dir_name}/{file_name}", "r") + + # ========================================================================== + # Stopping power for protons + # ========================================================================== + if "stopping_power" in file: + self.stopping_power = file["stopping_power"]["total_stopping_power"][()] + self.stopping_power_energy_grid = file["stopping_power"]["energy"][()] + elif simulation.settings.csda: + raise ValueError(f"CSDA cannot be used if no stopping power is provided for nuclide {self.name}") + + # Only CSDA data available - no nuclear rxn xs + if "proton_reactions" not in file: + # Zero out all xs arrays + xs_energy = np.array([0, 1.0e10]) + self.proton_xs_energy_grid = xs_energy + + self.proton_total_xs = np.zeros_like(self.proton_xs_energy_grid) + self.proton_elastic_xs = np.zeros_like(self.proton_xs_energy_grid) + self.proton_inelastic_xs = np.zeros_like(self.proton_xs_energy_grid) + + file.close() + return + + rx_names = [ + "elastic_scattering", + "inelastic_scattering", + "capture", + ] + + # The reaction MTs + MTs = {} + for name in rx_names: + if name not in file["proton_reactions"]: + MTs[name] = [] + continue + + MTs[name] = [ + x for x in file[f"proton_reactions/{name}"] if x.startswith("MT") + ] + + # ========================================================================== + # Reaction XS + # ========================================================================== + + # Energy grid + xs_energy = file["proton_reactions/xs_energy_grid"][()] * 1e6 # MeV to eV + self.proton_xs_energy_grid = xs_energy + + # The total XS + self.proton_total_xs = np.zeros_like(self.proton_xs_energy_grid) + self.proton_elastic_xs = np.zeros_like(self.proton_xs_energy_grid) + self.proton_inelastic_xs = np.zeros_like(self.proton_xs_energy_grid) + self.proton_capture_xs = np.zeros_like(self.proton_xs_energy_grid) + + xs_containers = [ + self.proton_elastic_xs, + self.proton_inelastic_xs, + self.proton_capture_xs, + ] + + for xs_container, rx_name in list(zip(xs_containers, rx_names)): + for MT in MTs[rx_name]: + xs = file[f"proton_reactions/{rx_name}/{MT}/xs"] + xs_container[xs.attrs["offset"] :] += xs[()] + + self.proton_total_xs = self.proton_elastic_xs + self.proton_inelastic_xs + self.proton_capture_xs + + # ========================================================================== + # The reactions + # ========================================================================== + + self.proton_elastic_scattering_reactions = [] + self.proton_inelastic_scattering_reactions = [] + self.proton_capture_reactions = [] + + rx_containers = [ + self.proton_elastic_scattering_reactions, + self.proton_inelastic_scattering_reactions, + self.proton_capture_reactions, + ] + rx_classes = [ + ProtonReactionElasticScattering, + ProtonReactionInelasticScattering, + ProtonReactionCapture, + ] + 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"proton_reactions/{rx_name}/{MT}"] + reaction = rx_class.from_h5_group(h5_group) + rx_container.append(reaction) + + file.close() + + ## TODO: UPDATE this to handle protons as well as neutrons def __repr__(self): text = "\n" text += f"Nuclide\n" diff --git a/mcdc/object_/proton_reaction.py b/mcdc/object_/proton_reaction.py new file mode 100644 index 000000000..57cda4f36 --- /dev/null +++ b/mcdc/object_/proton_reaction.py @@ -0,0 +1,392 @@ +from typing import Annotated +import numpy as np +from numpy import float64 +from numpy.typing import NDArray + +#### + +import mcdc.object_.distribution as distribution + +from mcdc.constant import ( + ANGLE_ISOTROPIC, + ANGLE_ENERGY_CORRELATED, + ANGLE_DISTRIBUTED, + INTERPOLATION_LINEAR, + INTERPOLATION_LOG, + PROTON_REACTION_ELASTIC_SCATTERING, + PROTON_REACTION_CAPTURE, + PROTON_REACTION_INELASTIC_SCATTERING, + REFERENCE_FRAME_COM, + REFERENCE_FRAME_LAB, + PARTICLE_NEUTRON, + PARTICLE_PROTON, +) +from mcdc.object_.base import ObjectPolymorphic +from mcdc.object_.distribution import ( + DistributionBase, + DistributionMultiTable, + DistributionLevelScattering, + DistributionEvaporation, + DistributionMaxwellian, + DistributionKalbachMann, + DistributionTabulatedEnergyAngle, + DistributionNBody, +) +from mcdc.object_.simulation import simulation +from mcdc.print_ import print_1d_array, print_error + +# ====================================================================================== +# Proton reaction base class +# ====================================================================================== + + +class ProtonReactionBase(ObjectPolymorphic): + # Annotations for Numba mode + label: str = "proton_reaction" + # + MT: int + xs: NDArray[float64] + xs_offset_: int # "xs_offset" ir reserved for "xs" + reference_frame: int + q_value: float64 + + def __init__(self, type_, MT, xs, xs_offset, reference_frame, q_value): + self.MT = MT + self.xs = xs + self.xs_offset_ = xs_offset + self.reference_frame = reference_frame + self.q_value = q_value + super().__init__(type_) + + def __repr__(self): + text = "\n" + text += f"{decode_type(self.type)}\n" + text += f" - ID: {self.ID}\n" + text += f" - MT: {self.MT}\n" + text += f" - XS {print_1d_array(self.xs)} barn\n" + text += f" - Reference frame: {decode_reference_frame(self.reference_frame)}\n" + text += f" - Q-value: {self.q_value}\n" + return text + + +def decode_type(type_): + if type_ == PROTON_REACTION_ELASTIC_SCATTERING: + return "Proton elastic scattering" + elif type_ == PROTON_REACTION_INELASTIC_SCATTERING: + return "Proton inelastic scattering" + elif type_ == PROTON_REACTION_CAPTURE: + return "Proton capture" + + +def decode_reference_frame(type_): + if type_ == REFERENCE_FRAME_LAB: + return "Laboratory" + elif type_ == REFERENCE_FRAME_COM: + return "Center of mass" + + +# ====================================================================================== +# Proton elastic scattering +# ====================================================================================== + + +class ProtonReactionElasticScattering(ProtonReactionBase): + # Annotations for Numba mode + label: str = "proton_elastic_scattering_reaction" + # + mu_table: DistributionBase + + def __init__(self, MT, xs, xs_offset, reference_frame, mu): + type_ = PROTON_REACTION_ELASTIC_SCATTERING + self.mu_table = mu + super().__init__(type_, MT, xs, xs_offset, reference_frame, 0.0) + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame, _ = set_basic_properties(h5_group) + _, mu = set_angular_distribution(h5_group["angular_cosine_distribution"]) + return cls(MT, xs, xs_offset, reference_frame, mu) + + def __repr__(self): + text = super().__repr__() + text += f" - Scattering cosine: {distribution.decode_type(self.mu_table.type)} [ID: {self.mu_table.ID}]\n" + return text + + +# ====================================================================================== +# Proton inelastic scattering +# ====================================================================================== + + +class ProtonReactionInelasticScattering(ProtonReactionBase): + # Annotations for Numba mode + label: str = "proton_inelastic_scattering_reaction" + # + multiplicity: int + angle_type: int + mu: DistributionBase + N_spectrum_probability_bin: int + N_spectrum: int + spectrum_probability_grid: NDArray[float64] + spectrum_probability: Annotated[ + NDArray[float64], ("N_spectrum_probability_bin", "N_spectrum") + ] + energy_spectra: list[DistributionBase] + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + q_value, + multiplicity, + angle_type, + mu, + spectrum_probability_grid, + spectrum_probability, + energy_spectra, + ): + type_ = PROTON_REACTION_INELASTIC_SCATTERING + super().__init__(type_, MT, xs, xs_offset, reference_frame, q_value) + + self.multiplicity = multiplicity + self.angle_type = angle_type + self.mu = mu + self.N_spectrum_probability_bin = len(spectrum_probability_grid) - 1 + self.N_spectrum = len(energy_spectra) + self.spectrum_probability_grid = spectrum_probability_grid + self.spectrum_probability = spectrum_probability + self.energy_spectra = energy_spectra + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame, q_value = set_basic_properties(h5_group) + multiplicity = int(h5_group["multiplicity"][()]) + + ang_type_str = h5_group["angular_cosine_distribution"].attrs.get("type", "isotropic") + if ang_type_str == "given_in_energy_distribution": + angle_type, mu = set_angular_distribution_from_kalbach_mann( + h5_group["energy_spectrum-1"] + ) + else: + angle_type, mu = set_angular_distribution( + h5_group["angular_cosine_distribution"] + ) + + spectrum_probability_grid = h5_group["spectrum_probability_grid"][()] * 1e6 + spectrum_probability = h5_group["spectrum_probability"][()] + energy_spectra = [ + set_energy_distribution(h5_group[name]) + for name in sorted(x for x in h5_group if x.startswith("energy_spectrum-")) + ] + + return cls(MT, xs, xs_offset, reference_frame, q_value, multiplicity, + angle_type, mu, spectrum_probability_grid, spectrum_probability, + energy_spectra) + + def __repr__(self): + text = super().__repr__() + if self.angle_type == ANGLE_ISOTROPIC: + text += f" - Scattering cosine: Isotropic\n" + elif self.angle_type == ANGLE_ENERGY_CORRELATED: + text += f" - Scattering cosine: Energy-correlated\n" + else: + text += f" - Scattering cosine: {distribution.decode_type(self.mu.type)} [ID: {self.mu.ID}]\n" + text += f" - Energy spectra\n" + text += f" - Probability energy grid {print_1d_array(self.spectrum_probability_grid)}\n" + for i in range(len(self.energy_spectra)): + text += f" - Spectrum {i+1}: {distribution.decode_type(self.energy_spectra[i])} [{print_1d_array(self.spectrum_probability[:,i])}] [ID: {self.energy_spectra[i].ID}]\n" + return text + + + +class ProtonReactionCapture(ProtonReactionBase): + # Annotations for Numba mode + label: str = "proton_capture_reaction" + + def __init__(self, MT, xs, xs_offset, reference_frame, q_value): + type_ = PROTON_REACTION_CAPTURE + super().__init__(type_, MT, xs, xs_offset, reference_frame, q_value) + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame, q_value = set_basic_properties(h5_group) + return cls(MT, xs, xs_offset, reference_frame, q_value) + + + +# ====================================================================================== +# Helper functions +# ====================================================================================== + + +def set_basic_properties(h5_group): + MT = int(h5_group.attrs["MT"][()]) + xs = h5_group["xs"][()] + xs_offset = h5_group["xs"].attrs["offset"] + reference_frame = h5_group["reference_frame"][()].decode("utf-8") + if reference_frame == "LAB": + reference_frame = REFERENCE_FRAME_LAB + elif reference_frame == "COM": + reference_frame = REFERENCE_FRAME_COM + q_value = h5_group["Q-value"][()] + return MT, xs, xs_offset, reference_frame, q_value + + +def set_angular_distribution(h5_group): + # Handle missing type attribute + if "type" not in h5_group.attrs: + mu_type = "isotropic" + else: + mu_type = h5_group.attrs["type"] + + if mu_type == "isotropic": + angle_type = ANGLE_ISOTROPIC + mu = simulation.distributions[0] + elif mu_type == "energy-correlated": + angle_type = ANGLE_ENERGY_CORRELATED + mu = simulation.distributions[0] + elif mu_type == "given_in_energy_distribution": + raise ValueError( + "set_angular_distribution called with given_in_energy_distribution; " + "use set_angular_distribution_from_kalbach_mann instead.") + elif mu_type == "tabulated": + angle_type = ANGLE_DISTRIBUTED + + # Check if data is in flattened format or subgroup format + if "energy" in h5_group: + # Flattened format + grid = h5_group[f"energy"][()] * 1e6 # MeV to eV + offset = h5_group[f"offset"][()] + value = h5_group[f"value"][()] + pdf = h5_group[f"pdf"][()] + else: + # Subgroup format: E_in_1, E_in_2, etc. + incident_energies = h5_group["incident_energies"][()] * 1e6 # MeV to eV + + # Collect all cosines and pdfs into flattened arrays + cosines_list = [] + pdf_list = [] + offset = np.zeros(len(incident_energies), dtype=np.int32) + + for i, energy in enumerate(incident_energies): + subgroup_name = f"E_in_{i + 1}" + if subgroup_name in h5_group: + subgroup = h5_group[subgroup_name] + if subgroup.attrs.get("type", "tabulated") == "tabulated": + cosines_list.extend(subgroup["cosines"][()]) + pdf_list.extend(subgroup["pdf"][()]) + else: + # Isotropic - use dummy values + cosines_list.extend([0.0]) # isotropic cosine + pdf_list.extend([1.0]) # uniform pdf + else: + # Missing subgroup - assume isotropic + cosines_list.extend([0.0]) + pdf_list.extend([1.0]) + + if i < len(incident_energies) - 1: + offset[i + 1] = len(cosines_list) + + grid = incident_energies + value = np.array(cosines_list) + pdf = np.array(pdf_list) + + mu = DistributionMultiTable(grid, offset, value, pdf) + + return angle_type, mu + +def set_angular_distribution_from_kalbach_mann(spectrum_group): + """ + Build a DistributionMultiTable for Kalbach-Mann angular sampling. + The 'value' array holds the angular slope 'a'. The transport kernel + uses these to sample cosines analytically via the Kalbach-Mann formula. + """ + grid = spectrum_group["energy"][()] * 1e6 # MeV to eV + offset = spectrum_group["offset"][()] + a = spectrum_group["angular_slope"][()] + pdf = spectrum_group["pdf"][()] + + mu = DistributionMultiTable(grid, offset, a, pdf) + return ANGLE_ENERGY_CORRELATED, mu + + +def set_energy_distribution(h5_group): + spectrum_type = h5_group.attrs["type"] + + if spectrum_type == "tabulated": + grid = h5_group[f"energy"][()] * 1e6 # MeV to eV + offset = h5_group[f"offset"][()] + value = h5_group[f"value"][()] * 1e6 # MeV to eV + pdf = h5_group[f"pdf"][()] / 1e6 # /MeV to /eV + energy_spectrum = DistributionMultiTable(grid, offset, value, pdf) + + elif spectrum_type == "level-scattering": + C1 = h5_group["C1"][()] * 1e6 # MeV to eV + C2 = h5_group["C2"][()] + + energy_spectrum = DistributionLevelScattering(C1, C2) + + elif spectrum_type == "evaporation": + energy = h5_group[f"temperature_energy_grid"][()] * 1e6 # MeV to eV + temperature = h5_group[f"temperature"][()] * 1e6 # MeV to eV + restriction_energy = h5_group[f"restriction_energy"][()] * 1e6 # MeV to eV + + energy_spectrum = DistributionEvaporation( + energy, temperature, restriction_energy + ) + + elif spectrum_type == "maxwellian": + energy = h5_group[f"temperature_energy_grid"][()] * 1e6 # MeV to eV + temperature = h5_group[f"temperature"][()] * 1e6 # MeV to eV + restriction_energy = h5_group[f"restriction_energy"][()] * 1e6 # MeV to eV + interpolation = h5_group[f"temperature_interpolation"][()].decode("utf-8") + if interpolation == "linear": + interpolation = INTERPOLATION_LINEAR + elif interpolation == "log": + interpolation = INTERPOLATION_LOG + + energy_spectrum = DistributionMaxwellian( + energy, temperature, restriction_energy, interpolation + ) + + elif spectrum_type == "kalbach-mann": + energy = h5_group[f"energy"][()] * 1e6 # MeV to eV + offset = h5_group[f"offset"][()] + + energy_out = h5_group[f"energy_out"][()] * 1e6 # MeV to eV + pdf = h5_group[f"pdf"][()] / 1e6 # /MeV to /eV + + precompound_factor = h5_group[f"precompound_factor"][()] + angular_slope = h5_group[f"angular_slope"][()] + + energy_spectrum = DistributionKalbachMann( + energy, offset, energy_out, pdf, precompound_factor, angular_slope + ) + + elif spectrum_type == "energy-angle-tabulated": + energy = h5_group[f"energy"][()] * 1e6 # MeV to eV + offset = h5_group[f"offset"][()] + + energy_out = h5_group[f"energy_out"][()] * 1e6 # MeV to eV + pdf = h5_group[f"pdf"][()] / 1e6 # /MeV to /eV + cosine_offset = h5_group[f"cosine_offset"][()] + + cosine = h5_group[f"cosine"][()] + cosine_pdf = h5_group[f"cosine_pdf"][()] + + energy_spectrum = DistributionTabulatedEnergyAngle( + energy, offset, energy_out, pdf, cosine_offset, cosine, cosine_pdf + ) + + elif spectrum_type == "N-body": + value = h5_group["value"][()] * 1e6 # MeV to eV + pdf = h5_group["pdf"][()] / 1e6 # /MeV to /eV + + energy_spectrum = DistributionNBody(value, pdf) + + else: + print_error(f"Unsupported energy spectrum of type {spectrum_type}") + + return energy_spectrum diff --git a/mcdc/object_/settings.py b/mcdc/object_/settings.py index 60a28d896..2caac29c9 100644 --- a/mcdc/object_/settings.py +++ b/mcdc/object_/settings.py @@ -44,6 +44,8 @@ class Settings(ObjectSingleton): time_boundary: float = np.inf output_name: str = "output" use_progress_bar: bool = True + csda: bool = False + csda_max_fractional_e_loss: float = 0.01 # Time census N_census: int = 1 diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py index 48e7e3433..6787de492 100644 --- a/mcdc/object_/simulation.py +++ b/mcdc/object_/simulation.py @@ -16,6 +16,7 @@ from mcdc.object_.material import MaterialBase from mcdc.object_.nuclide import Nuclide from mcdc.object_.neutron_reaction import NeutronReactionBase + from mcdc.object_.proton_reaction import ProtonReactionBase from mcdc.object_.source import Source from mcdc.object_.surface import Surface from mcdc.object_.tally import Tally @@ -64,6 +65,7 @@ class Simulation(ObjectSingleton): nuclides: list[Nuclide] neutron_reactions: list[NeutronReactionBase] sources: list[Source] + proton_reactions: list[ProtonReactionBase] # Geometry cells: list[Cell] @@ -150,6 +152,7 @@ def __init__(self): self.nuclides = [] self.neutron_reactions = [] self.sources = [] + self.proton_reactions = [] # Geometry self.cells = [] diff --git a/mcdc/transport/physics/__init__.py b/mcdc/transport/physics/__init__.py index 66133009c..c0e782318 100644 --- a/mcdc/transport/physics/__init__.py +++ b/mcdc/transport/physics/__init__.py @@ -4,6 +4,9 @@ neutron_production_xs, collision_distance, collision, + csda_distance, + csda_edep, ) import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron +import mcdc.transport.physics.proton as proton diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py index ef3cef349..29ace02a6 100644 --- a/mcdc/transport/physics/interface.py +++ b/mcdc/transport/physics/interface.py @@ -1,4 +1,5 @@ import math +import numpy as np from numba import njit @@ -7,6 +8,9 @@ import mcdc.transport.rng as rng import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron +import mcdc.transport.physics.proton as proton + +import mcdc.mcdc_get as mcdc_get from mcdc.constant import * @@ -22,6 +26,8 @@ def particle_speed(particle_container, simulation, data): return neutron.particle_speed(particle_container, simulation, data) elif particle["particle_type"] == PARTICLE_ELECTRON: return electron.particle_speed(particle_container, simulation, data) + elif particle["particle_type"] == PARTICLE_PROTON: + return proton.particle_speed(particle_container, simulation, data) return -1.0 @@ -37,6 +43,8 @@ def macro_xs(reaction_type, particle_container, simulation, data): return neutron.macro_xs(reaction_type, particle_container, simulation, data) elif particle["particle_type"] == PARTICLE_ELECTRON: return electron.macro_xs(reaction_type, particle_container, simulation, data) + elif particle["particle_type"] == PARTICLE_PROTON: + return proton.macro_xs(reaction_type, particle_container, simulation, data) return -1.0 @@ -65,6 +73,8 @@ def collision_distance(particle_container, simulation, data): SigmaT = macro_xs(NEUTRON_REACTION_TOTAL, particle_container, simulation, data) elif particle["particle_type"] == PARTICLE_ELECTRON: SigmaT = macro_xs(ELECTRON_REACTION_TOTAL, particle_container, simulation, data) + elif particle["particle_type"] == PARTICLE_PROTON: + SigmaT = macro_xs(PROTON_REACTION_TOTAL, particle_container, simulation, data) # Vacuum material? if SigmaT == 0.0: @@ -84,3 +94,57 @@ def collision(particle_container, collision_data_container, program, data): neutron.collision(particle_container, collision_data_container, program, data) elif particle["particle_type"] == PARTICLE_ELECTRON: electron.collision(particle_container, collision_data_container, program, data) + elif particle["particle_type"] == PARTICLE_PROTON: + proton.collision(particle_container, collision_data_container, program, data) + + +# ====================================================================================== +# Continuous Slowing Down Approximation +# ====================================================================================== + + +@njit +def csda_distance(particle_container, simulation, data): + particle = particle_container[0] + material = simulation["native_materials"][particle["material_ID"]] + E = particle["E"] + total_rho = 0.0 + total_dedx = 0.0 + + for i in range(material["N_nuclide"]): + nuclide_ID = int(mcdc_get.native_material.nuclide_IDs(i, material, data)) + nuclide = simulation["nuclides"][nuclide_ID] + + if not material["stopping_power_provided"]: + dedx_values = mcdc_get.nuclide.stopping_power_all(nuclide, data) + dedx_energies = mcdc_get.nuclide.stopping_power_energy_grid_all(nuclide, data) + dedx = np.interp(E / 1e6, dedx_energies, dedx_values) + total_dedx += dedx * 1e6 + + atomic_mass = nuclide["atomic_weight_ratio"] + nuclide_density = mcdc_get.native_material.nuclide_densities(i, material, data) + density_gcm3 = nuclide_density * 1e24 * atomic_mass / (6.022e23) + total_rho += density_gcm3 + + if material["stopping_power_provided"]: + dedx_values = mcdc_get.native_material.stopping_power_all(material, data) + dedx_energies = mcdc_get.native_material.stopping_power_energy_grid_all(material, data) + dedx = np.interp(E / 1e6, dedx_energies, dedx_values) + total_dedx = dedx * 1e6 + + + max_fractional_e_loss = simulation["settings"]["csda_max_fractional_e_loss"] + return max_fractional_e_loss * E / total_dedx / total_rho + + +@njit +def csda_edep(particle_container, collision_data_container, distance, simulation, data): + particle = particle_container[0] + if particle["particle_type"] == PARTICLE_NEUTRON: + raise ValueError("CSDA not supported for neutrons") + if particle["particle_type"] == PARTICLE_ELECTRON: + raise ValueError("CSDA not supported for electrons") + if particle["particle_type"] == PARTICLE_PROTON: + proton.csda_edep( + particle_container, collision_data_container, distance, simulation, data + ) diff --git a/mcdc/transport/physics/proton/__init__.py b/mcdc/transport/physics/proton/__init__.py new file mode 100644 index 000000000..6fed2a1e0 --- /dev/null +++ b/mcdc/transport/physics/proton/__init__.py @@ -0,0 +1,8 @@ +from .interface import ( + particle_speed, + macro_xs, + collision, + csda_edep, +) +import mcdc.transport.physics.proton.native as native +import mcdc.transport.physics.proton.multigroup as multigroup diff --git a/mcdc/transport/physics/proton/interface.py b/mcdc/transport/physics/proton/interface.py new file mode 100644 index 000000000..34b81c01f --- /dev/null +++ b/mcdc/transport/physics/proton/interface.py @@ -0,0 +1,49 @@ +from numba import njit + +#### + +import mcdc.transport.physics.proton.multigroup as multigroup +import mcdc.transport.physics.proton.native as native +import mcdc.transport.util as util + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container, simulation, data): + return native.particle_speed(particle_container) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, simulation, data): + return native.macro_xs(reaction_type, particle_container, simulation, data) + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, collision_data_container, program, data): + simulation = util.access_simulation(program) + native.collision(particle_container, collision_data_container, program, data) + + +# ====================================================================================== +# Continuous Slowing Down Approximation +# ====================================================================================== + + +@njit +def csda_edep(particle_container, collision_data_container, distance, simulation, data): + native.csda_edep( + particle_container, collision_data_container, distance, simulation, data + ) diff --git a/mcdc/transport/physics/proton/multigroup.py b/mcdc/transport/physics/proton/multigroup.py new file mode 100644 index 000000000..b57c5b3b6 --- /dev/null +++ b/mcdc/transport/physics/proton/multigroup.py @@ -0,0 +1,170 @@ +import numpy as np +import math + +from numba import njit + +#### + +import mcdc.mcdc_get as mcdc_get +import mcdc.numba_types as type_ +import mcdc.transport.particle as particle_module +import mcdc.transport.particle_bank as particle_bank_module +import mcdc.transport.rng as rng +import mcdc.transport.util as util + +from mcdc.constant import ( + PI, + PROTON_REACTION_TOTAL, + PROTON_REACTION_ELASTIC_SCATTERING, + PROTON_REACTION_INELASTIC_SCATTERING, + PROTON_REACTION_CAPTURE, +) +from mcdc.transport.physics.util import scatter_direction +from mcdc.transport.distribution import sample_isotropic_direction + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container, simulation, data): + particle = particle_container[0] + material = simulation["multigroup_materials"][particle["material_ID"]] + return mcdc_get.multigroup_material.mgxs_speed(particle["g"], material, data) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, simulation, data): + particle = particle_container[0] + material = simulation["multigroup_materials"][particle["material_ID"]] + g = particle["g"] + + if reaction_type == PROTON_REACTION_TOTAL: + return mcdc_get.multigroup_material.mgxs_total(g, material, data) + elif reaction_type == PROTON_REACTION_ELASTIC_SCATTERING: + return mcdc_get.multigroup_material.mgxs_scatter(g, material, data) + return 0.0 + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, collision_data_container, program, data): + simulation = util.access_simulation(program) + particle = particle_container[0] + + # Get the reaction cross-sections + SigmaT = macro_xs(PROTON_REACTION_TOTAL, particle_container, simulation, data) + SigmaS = macro_xs( + PROTON_REACTION_ELASTIC_SCATTERING, particle_container, simulation, data + ) + SigmaC = macro_xs(PROTON_REACTION_CAPTURE, particle_container, simulation, data) + + # Implicit capture + if simulation["implicit_capture"]["active"]: + particle["w"] *= (SigmaT - SigmaC) / SigmaT + SigmaT -= SigmaC + + # Sample reaction type and perform the reaction + xi = rng.lcg(particle_container) * SigmaT + total = SigmaS + if total > xi: + scattering(particle_container, program, data) + else: + particle["alive"] = False + + +# ====================================================================================== +# Reactions +# ====================================================================================== + + +@njit +def scattering(particle_container, program, data): + simulation = util.access_simulation(program) + + # Particle attributes + particle = particle_container[0] + g = particle["g"] + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + + # Material attributes + material = simulation["multigroup_materials"][particle["material_ID"]] + G = material["G"] + + # Kill the current particle + particle["alive"] = False + + # Adjust production and product weights if weighted emission + weight_production = 1.0 + weight_product = particle["w"] + if simulation["weighted_emission"]["active"]: + weight_target = simulation["weighted_emission"]["weight_target"] + weight_production = particle["w"] / weight_target + weight_product = weight_target + + + # TODO: make this better for protons, add secondary particle generation to non-MG materials + # Get number of secondaries + nu_s = mcdc_get.multigroup_material.mgxs_nu_s(g, material, data) + N = int(math.floor(weight_production * nu_s + rng.lcg(particle_container))) + + # Set up secondary partice container + particle_container_new = util.local_array(1, type_.particle_data) + particle_new = particle_container_new[0] + + # Create the secondaries + for n in range(N): + # Set default attributes + particle_module.copy_as_child(particle_container_new, particle_container) + + # Set weight + particle_new["w"] = weight_product + + # Sample scattering angle + mu0 = 2.0 * rng.lcg(particle_container_new) - 1.0 + + # Scatter direction + azi = 2.0 * PI * rng.lcg(particle_container_new) + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) + particle_new["ux"] = ux_new + particle_new["uy"] = uy_new + particle_new["uz"] = uz_new + + # Get outgoing spectrum + stride = material["G"] + start = material["mgxs_chi_s_offset"] + g * stride + chi_s = data[start : start + stride] + # Above is equivalent to: chi_s = mcdc_get.multigroup_material.mgxs_chi_s_vector(g, material, data) + + # Sample outgoing energy + xi = rng.lcg(particle_container_new) + total = 0.0 + for g_out in range(G): + total += chi_s[g_out] + if total > xi: + break + particle_new["g"] = g_out + + # Bank, but keep it if it is the last particle + if n == N - 1: + particle["alive"] = True + particle["ux"] = particle_new["ux"] + particle["uy"] = particle_new["uy"] + particle["uz"] = particle_new["uz"] + particle["g"] = particle_new["g"] + particle["E"] = particle_new["E"] + particle["w"] = particle_new["w"] + else: + particle_bank_module.bank_active_particle(particle_container_new, program) diff --git a/mcdc/transport/physics/proton/native.py b/mcdc/transport/physics/proton/native.py new file mode 100644 index 000000000..50e449e58 --- /dev/null +++ b/mcdc/transport/physics/proton/native.py @@ -0,0 +1,839 @@ +import math +import numpy as np +from numba import njit +import time + +#### + +import mcdc.mcdc_get as mcdc_get +import mcdc.numba_types as type_ +import mcdc.transport.particle as particle_module +import mcdc.transport.particle_bank as particle_bank_module +import mcdc.transport.rng as rng +import mcdc.transport.util as util + +from mcdc.constant import ( + ANGLE_DISTRIBUTED, + ANGLE_ENERGY_CORRELATED, + ANGLE_ISOTROPIC, + BOLTZMANN_K, + THERMAL_THRESHOLD_FACTOR, + LIGHT_SPEED, + PROTON_MASS, + PI, + PI_HALF, + PI_SQRT, + PROTON_REACTION_TOTAL, + PROTON_REACTION_ELASTIC_SCATTERING, + PROTON_REACTION_CAPTURE, + PROTON_REACTION_INELASTIC_SCATTERING, + REFERENCE_FRAME_COM, + PARTICLE_ELECTRON, + PARTICLE_NEUTRON, + PARTICLE_PROTON, + PROTON_CUTOFF_ENERGY, +) +from mcdc.transport.data import evaluate_data +from mcdc.transport.distribution import ( + sample_correlated_distribution_with_scale, + sample_distribution_with_scale, + sample_isotropic_cosine, + sample_isotropic_direction, + sample_multi_table, + sample_kalbach_mann, +) +from mcdc.transport.physics.util import ( + evaluate_proton_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 = PROTON_MASS + return LIGHT_SPEED * math.sqrt(E * (E + 2.0 * mass)) / (E + mass) + + +@njit +def particle_energy_from_speed(speed): + beta = speed / LIGHT_SPEED + gamma = 1.0 / math.sqrt(1.0 - beta * beta) + mass = PROTON_MASS + return mass * (gamma - 1.0) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, simulation, data): + particle = particle_container[0] + material = simulation["native_materials"][particle["material_ID"]] + E = particle["E"] + + total = 0.0 + + for i in range(material["N_nuclide"]): + nuclide_ID = int(mcdc_get.native_material.nuclide_IDs(i, material, data)) + nuclide = simulation["nuclides"][nuclide_ID] + + nuclide_density = mcdc_get.native_material.nuclide_densities(i, material, data) + xs = total_micro_xs(reaction_type, E, nuclide, data) + + total += nuclide_density * xs + + return total + + +@njit +def total_micro_xs(reaction_type, E, nuclide, data): + + idx, E0, E1 = evaluate_proton_xs_energy_grid(E, nuclide, data) + if reaction_type == PROTON_REACTION_TOTAL: + xs0 = mcdc_get.nuclide.proton_total_xs(idx, nuclide, data) + xs1 = mcdc_get.nuclide.proton_total_xs(idx + 1, nuclide, data) + + elif reaction_type == PROTON_REACTION_ELASTIC_SCATTERING: + xs0 = mcdc_get.nuclide.proton_elastic_xs(idx, nuclide, data) + xs1 = mcdc_get.nuclide.proton_elastic_xs(idx + 1, nuclide, data) + + # total_elastic_scattering_xs = 0 + # # Get total elastic scattering xs for all possible elastic scattering rxns + # for i in range(nuclide["N_proton_elastic_scattering_reaction"]): + # reaction_ID = int( + # mcdc_get.nuclide.proton_elastic_scattering_reaction_IDs(i, nuclide, data) + # ) + # reaction = simulation["proton_elastic_scattering_reactions"][reaction_ID] + # reaction_base_ID = reaction["parent_ID"] + # reaction_base = simulation["proton_reactions"][reaction_base_ID] + # xs = reaction_micro_xs(E, reaction_base, nuclide, data) + # total_elastic_scattering_xs += xs + + # return total_elastic_scattering_xs + + elif reaction_type == PROTON_REACTION_INELASTIC_SCATTERING: + xs0 = mcdc_get.nuclide.proton_inelastic_xs(idx, nuclide, data) + xs1 = mcdc_get.nuclide.proton_inelastic_xs(idx + 1, nuclide, data) + + # print(f'all inelastic xs = {mcdc_get.nuclide.proton_inelastic_xs_all(nuclide, data)}') + # print(f'idx = {idx}, xs0 = {xs0}, xs1 = {xs1}') + # raise ValueError("stop") + + # total_inelastic_scattering_xs = 0 + # # Get total inelastic scattering xs for all possible inelastic scattering rxns + # for i in range(nuclide["N_proton_inelastic_scattering_reaction"]): + # reaction_ID = int( + # mcdc_get.nuclide.proton_inelastic_scattering_reaction_IDs(i, nuclide, data) + # ) + # reaction = simulation["proton_inelastic_scattering_reactions"][reaction_ID] + # reaction_base_ID = reaction["parent_ID"] + # reaction_base = simulation["proton_reactions"][reaction_base_ID] + # xs = reaction_micro_xs(E, reaction_base, nuclide, data) + # total_inelastic_scattering_xs += xs + + # return total_inelastic_scattering_xs + elif reaction_type == PROTON_REACTION_CAPTURE: + xs0 = mcdc_get.nuclide.proton_capture_xs(idx, nuclide, data) + xs1 = mcdc_get.nuclide.proton_capture_xs(idx + 1, nuclide, data) + # total_capture_xs = 0 + # # Get total capture xs for all possible capture rxns + # for i in range(nuclide["N_proton_capture_reaction"]): + # reaction_ID = int( + # mcdc_get.nuclide.proton_capture_reaction_IDs(i, nuclide, data) + # ) + # reaction = simulation["proton_capture_reactions"][reaction_ID] + # reaction_base_ID = reaction["parent_ID"] + # reaction_base = simulation["proton_reactions"][reaction_base_ID] + # xs = reaction_micro_xs(E, reaction_base, nuclide, data) + # total_capture_xs += xs + + # return total_capture_xs + + else: + # Should be unreachable + xs0 = 0.0 + xs1 = 0.0 + return linear_interpolation(E, E0, E1, xs0, xs1) + + +@njit +def reaction_micro_xs(E, reaction_base, nuclide, data): + idx, E0, E1 = evaluate_proton_xs_energy_grid(E, nuclide, data) + + # Apply offset + offset = reaction_base["xs_offset_"] + if idx < offset: + return 0.0 + else: + idx -= offset + + xs0 = mcdc_get.proton_reaction.xs(idx, reaction_base, data) + xs1 = mcdc_get.proton_reaction.xs(idx + 1, reaction_base, data) + return linear_interpolation(E, E0, E1, xs0, xs1) + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, collision_data_container, program, data): + simulation = util.access_simulation(program) + particle = particle_container[0] + collision_data = collision_data_container[0] + material = simulation["native_materials"][particle["material_ID"]] + + # Particle properties + E = particle["E"] + + # Check for cutoff energy + if E <= PROTON_CUTOFF_ENERGY: + collision_data["energy_deposition"] += E * particle["w"] + particle["alive"] = False + particle["E"] = 0.0 + return + + # ================================================================================== + # Sample colliding nuclide + # ================================================================================== + + SigmaT = macro_xs(PROTON_REACTION_TOTAL, particle_container, simulation, data) + + # TODO: add implicit capture for protons + xi = rng.lcg(particle_container) * SigmaT + total = 0.0 + for i in range(material["N_nuclide"]): + nuclide_ID = int(mcdc_get.native_material.nuclide_IDs(i, material, data)) + nuclide = simulation["nuclides"][nuclide_ID] + + nuclide_density = mcdc_get.native_material.nuclide_densities(i, material, data) + sigmaT = total_micro_xs(PROTON_REACTION_TOTAL, E, nuclide, data) + + SigmaT_nuclide = nuclide_density * sigmaT + total += SigmaT_nuclide + + if total > xi: + break + + + # ================================================================================== + # Sample and perform reaction + # ================================================================================== + + sigma_elastic = total_micro_xs(PROTON_REACTION_ELASTIC_SCATTERING, E, nuclide, data) + sigma_inelastic = total_micro_xs(PROTON_REACTION_INELASTIC_SCATTERING, E, nuclide, data) + sigma_capture = total_micro_xs(PROTON_REACTION_CAPTURE, E, nuclide, data) + xi = rng.lcg(particle_container) * sigmaT + + # Elastic scattering + total = sigma_elastic + if xi < total: + # Sample the actual reaction from the group + total -= sigma_elastic + for i in range(nuclide["N_proton_elastic_scattering_reaction"]): + reaction_ID = int( + mcdc_get.nuclide.proton_elastic_scattering_reaction_IDs( + i, nuclide, data + ) + ) + reaction = simulation["proton_elastic_scattering_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = simulation["proton_reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, nuclide, data) + + # Execute the reaction + if xi < total: + elastic_scattering( + reaction, + particle_container, + collision_data_container, + nuclide, + simulation, + data, + ) + return + + # Capture + if not simulation["implicit_capture"]["active"]: + # print(f'particle being captured') + sigma_capture = total_micro_xs(PROTON_REACTION_CAPTURE, E, nuclide, data) + total += sigma_capture + if xi < total: + # Sample the actual reaction from the group + total -= sigma_capture + for i in range(nuclide["N_proton_capture_reaction"]): + reaction_ID = int(mcdc_get.nuclide.proton_capture_reaction_IDs(i, nuclide, data)) + reaction = simulation["proton_capture_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = simulation["proton_reactions"][reaction_base_ID] + xs = reaction_micro_xs(E, reaction_base, nuclide, data) + total += xs + + # Execute the reaction + if xi < total: + capture( + reaction, + particle_container, + collision_data_container, + nuclide, + simulation, + data + ) + + + # Inelastic scattering + total += sigma_inelastic + if xi < total: + # Sample the actual reaction from the group + total -= sigma_inelastic + for i in range(nuclide["N_proton_inelastic_scattering_reaction"]): + reaction_ID = int( + mcdc_get.nuclide.proton_inelastic_scattering_reaction_IDs(i, nuclide, data) + ) + reaction = simulation["proton_inelastic_scattering_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = simulation["proton_reactions"][reaction_base_ID] + xs = reaction_micro_xs(E, reaction_base, nuclide, data) + total += xs + + # Execute the reaction + if xi < total: + inelastic_scattering( + reaction, + particle_container, + collision_data_container, + nuclide, + program, + data, + ) + return + + +# ====================================================================================== +# Continous Slowing Down Approximation (CSDA) +# ====================================================================================== + + +@njit +def csda_edep(particle_container, collision_data_container, distance, simulation, data): + particle = particle_container[0] + collision_data = collision_data_container[0] + material = simulation["native_materials"][particle["material_ID"]] + E = particle["E"] + + # Check for cutoff energy + if E <= PROTON_CUTOFF_ENERGY: + collision_data["energy_deposition"] += E * particle["w"] + particle["alive"] = False + particle["E"] = 0.0 + return + + average_A, average_Z, total_stopping_power, total_rho_gcm3 = calculate_total_stopping_power(particle_container, simulation, data) + energy_loss = total_stopping_power * total_rho_gcm3 * distance + + # Range straggling - modify energy loss to have some slight variations + # TODO: Insert different thickness regimes to sample from (e.g. Bohr, Landau, Vavilov) + # TODO: Make this part use rng state instead of np.random.normal? + energy_straggling_variance = 0.1569 * total_rho_gcm3 * average_Z / average_A * distance + energy_straggling_modifier = np.random.normal(loc=0.0, scale=np.sqrt(energy_straggling_variance)) + energy_loss += energy_straggling_modifier + particle["E"] -= energy_loss + collision_data["energy_deposition"] += energy_loss * particle["w"] + + if energy_loss * particle["w"] <= 0.0: + print(f'total density = {total_rho_gcm3}') + print(f'stopping_power = {total_stopping_power}') + print(f'distance = {distance}') + print(f'energy_loss = {energy_loss * particle["w"]}') + raise ValueError('negative energy loss') + + X0 = material["radiation_length"] + + # Angular scattering according to MCS theory + phi, theta = sample_mcs_angle(particle["E"], distance, total_rho_gcm3, X0) + + rotate_direction(particle, phi, theta) + + return + + +# ====================================================================================== +# Capture +# ====================================================================================== + + +# TODO: add secondaries from capture rxns +@njit +def capture( + reaction, particle_container, collision_data_container, nuclide, simulation, data +): + particle = particle_container[0] + collision_data = collision_data_container[0] + + reaction_base_ID = reaction["parent_ID"] + reaction_base = simulation["proton_reactions"][reaction_base_ID] + + # Terminate the particle + particle["alive"] = False + + # Energy deposition + E = particle["E"] + q_value = reaction_base["q_value"] * 1e6 + collision_data["energy_deposition"] += (E + q_value) * particle["w"] + + + +# ====================================================================================== +# Elastic scattering +# ====================================================================================== + + +@njit +def elastic_scattering( + reaction, particle_container, collision_data_container, nuclide, simulation, data +): + + # print(f'reaction = {repr(reaction)}') + # print(f'{reaction.dtype.names}') + # print(f'{reaction["mu_table_ID"]}, {reaction["ID"]}, {reaction["parent_ID"]}') + + # print(f'particle undergoing elastic scattering') + particle = particle_container[0] + collision_data = collision_data_container[0] + + # Particle attributes + E = particle["E"] + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + + # Energy deposition + collision_data["energy_deposition"] += E * particle["w"] + + # Note: Q-value is zero in elastic scattering + + # Sample nucleus thermal velocity + A = nuclide["atomic_weight_ratio"] + temperature = nuclide["temperature"] + if E > THERMAL_THRESHOLD_FACTOR * BOLTZMANN_K * temperature: + Vx = 0.0 + Vy = 0.0 + Vz = 0.0 + else: + Vx, Vy, Vz = sample_nucleus_velocity(A, particle_container) + + # ========================================================================= + # COM kinematics + # ========================================================================= + + # Particle speed + speed = particle_speed(particle_container) + + # Proton velocity - LAB + vx = speed * ux + vy = speed * uy + vz = speed * uz + + # COM velocity + COM_x = (vx + A * Vx) / (1.0 + A) + COM_y = (vy + A * Vy) / (1.0 + A) + COM_z = (vz + A * Vz) / (1.0 + A) + + # Proton velocity - COM + vx = vx - COM_x + vy = vy - COM_y + vz = vz - COM_z + + # Proton speed - COM + speed = math.sqrt(vx * vx + vy * vy + vz * vz) + + # Proton initial direction - COM + ux = vx / speed + uy = vy / speed + uz = vz / speed + + # # Sample the scattering cosine from the multi-PDF distribution + # print(f'simulation = {simulation}, names = {simulation.dtype.names}') + # print(f'reaction = {reaction}, names = {reaction.dtype.names}') + multi_table = simulation["multi_table_distributions"][reaction["mu_table_ID"]] + mu0 = sample_multi_table(E, particle_container, multi_table, simulation, data) + + # Scatter the direction in COM + azi = 2.0 * PI * rng.lcg(particle_container) + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) + + # Proton final velocity - COM + vx = speed * ux_new + vy = speed * uy_new + vz = speed * uz_new + + # ========================================================================= + # COM to LAB + # ========================================================================= + + # Final velocity - LAB + vx = vx + COM_x + vy = vy + COM_y + vz = vz + COM_z + + # Final energy - LAB + speed = math.sqrt(vx * vx + vy * vy + vz * vz) + particle["E"] = particle_energy_from_speed(speed) + + # Final direction - LAB + particle["ux"] = vx / speed + particle["uy"] = vy / speed + particle["uz"] = vz / speed + + # Subtract outgoing energy from energy deposition + collision_data["energy_deposition"] -= particle["E"] * particle["w"] + + +@njit +def sample_nucleus_velocity(A, particle_container): + particle = particle_container[0] + + # Particle speed + speed = particle_speed(particle_container) + + # Maxwellian parameter + beta = math.sqrt(2.0659834e-11 * A) + # The constant above is + # (1.674927471e-27 kg) / (1.38064852e-19 cm^2 kg s^-2 K^-1) / (293.6 K)/2 + + # Sample nuclide speed candidate V_tilda and + # nuclide-proton polar cosine candidate mu_tilda via + # rejection sampling + y = beta * speed + while True: + if rng.lcg(particle_container) < 2.0 / (2.0 + PI_SQRT * y): + x = math.sqrt( + -math.log(rng.lcg(particle_container) * rng.lcg(particle_container)) + ) + else: + cos_val = math.cos(PI_HALF * rng.lcg(particle_container)) + x = math.sqrt( + -math.log(rng.lcg(particle_container)) + - math.log(rng.lcg(particle_container)) * cos_val * cos_val + ) + V_tilda = x / beta + mu_tilda = 2.0 * rng.lcg(particle_container) - 1.0 + + # Accept candidate V_tilda and mu_tilda? + if rng.lcg(particle_container) > math.sqrt( + speed * speed + V_tilda * V_tilda - 2.0 * speed * V_tilda * mu_tilda + ) / (speed + V_tilda): + break + + # Set nuclide velocity - LAB + azi = 2.0 * PI * rng.lcg(particle_container) + ux, uy, uz = scatter_direction( + particle["ux"], particle["uy"], particle["uz"], mu_tilda, azi + ) + Vx = ux * V_tilda + Vy = uy * V_tilda + Vz = uz * V_tilda + + return Vx, Vy, Vz + + +# ====================================================================================== +# Inelastic scattering +# ====================================================================================== + + +# TODO: make inelastic scattering actually produce secondaries +@njit +def inelastic_scattering( + reaction, particle_container, collision_data_container, nuclide, program, data +): + # """ + # Proton intelastic scattering with secondary particle production. + + # Samples: + # 1. Outgoing proton from proton_reactions/inelastic_scattering/MT-005 + # 2. Secondary particles from secondary_particles/ZAP_x/MT-005 + # """ + # print(f'particle undergoing inelastic scattering') + + simulation = util.access_simulation(program) + particle = particle_container[0] + collision_data = collision_data_container[0] + + reaction_base_ID = reaction["parent_ID"] + reaction_base = simulation["proton_reactions"][reaction_base_ID] + + # Particle attributes + E = particle["E"] + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + w = particle["w"] + + # Kill the incident proton + particle["alive"] = False + + # Q-value energy available + q_value = reaction_base["q_value"] * 1e6 + total_energy = E + q_value + + # =========================================================================== + # Sample outgoing proton + # =========================================================================== + + # Number of outgoing protons and spectra + N_proton = reaction["multiplicity"] + N_spectrum = reaction["N_spectrum"] + use_all_spectrum = N_proton == N_spectrum + + # Set up secondary particle container + particle_container_new = util.local_array(1, type_.particle_data) + particle_new = particle_container_new[0] + + # Energy deposition (will be adjusted as we create secondaries) + collision_data["energy_deposition"] += total_energy * w + + # Create outgoing protons + for n in range(N_proton): + # Set default attributes (copy incident proton) + particle_module.copy_as_child(particle_container_new, particle_container) + + + + # ============================================================================== + # Sample angle (if not energy-correlated) + # ============================================================================== + + angle_type = reaction["angle_type"] + if angle_type == ANGLE_ENERGY_CORRELATED: + pass + elif angle_type == ANGLE_ISOTROPIC: + mu = sample_isotropic_cosine(particle_container_new) + elif angle_type == ANGLE_DISTRIBUTED: + distribution_base = simulation["distributions"][reaction["mu_ID"]] + multi_table = simulation["multi_table_distributions"][ + distribution_base["child_ID"] + ] + + mu = sample_multi_table(E, particle_container, multi_table, simulation, data) + + # ============================================================================== + # Sample energy (also angle if correlated) + # ============================================================================== + + # Get energy spectrum + if use_all_spectrum: + ID = int( + mcdc_get.proton_inelastic_scattering_reaction.energy_spectrum_IDs( + n, reaction, data + ) + ) + spectrum_base = simulation["distributions"][ID] + else: + offset = reaction["spectrum_probability_grid_offset"] + length = reaction["spectrum_probability_grid_length"] + probability_grid = data[offset : offset + length] + probability_idx = find_bin(E, probability_grid) + xi = rng.lcg(particle_container_new) + total = 0.0 + for j in range(N_spectrum): + probability = mcdc_get.proton_inelastic_scattering_reaction.spectrum_probability( + probability_idx, j, reaction, data + ) + total += probability + if xi < total: + ID = int( + mcdc_get.proton_inelastic_scattering_reaction.energy_spectrum_IDs( + j, reaction, data + ) + ) + spectrum_base = simulation["distributions"][ID] + break + + # Sample energy + if not angle_type == ANGLE_ENERGY_CORRELATED: + E_new = sample_distribution_with_scale( + E, spectrum_base, particle_container_new, simulation, data + ) + else: + E_new, mu = sample_correlated_distribution_with_scale( + E, spectrum_base, particle_container_new, simulation, data + ) + + # ============================================================================== + # Frame transformation + # ============================================================================== + + reference_frame = reaction_base["reference_frame"] + if reference_frame == REFERENCE_FRAME_COM: + A = nuclide["atomic_weight_ratio"] + mu_COM = mu + E_COM = E_new + + E_new = ( + E_COM + (E + 2 * mu_COM * (A + 1) * math.sqrt(E * E_COM)) / (A + 1) ** 2 + ) + mu = mu_COM * math.sqrt(E_COM / E_new) + math.sqrt(E / E_new) / (A + 1) + + azi = 2.0 * PI * rng.lcg(particle_container_new) + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu, azi) + + # Now the secondary angle and energy are finalized + particle_new["ux"] = ux_new + particle_new["uy"] = uy_new + particle_new["uz"] = uz_new + particle_new["E"] = E_new + particle_new["particle_type"] = PARTICLE_PROTON + + # Subtract outgoing energy from energy deposition + collision_data["energy_deposition"] -= particle_new["E"] * particle_new["w"] + + # ============================================================================== + # Bank the new particle + # ============================================================================== + + # Keep it if it is the last particle + if n == N_proton - 1: + particle["alive"] = True + particle["ux"] = particle_new["ux"] + particle["uy"] = particle_new["uy"] + particle["uz"] = particle_new["uz"] + particle["E"] = particle_new["E"] + particle["particle_type"] = PARTICLE_PROTON + else: + particle_bank_module.bank_active_particle(particle_container_new, program) + + # =========================================================================== + # 2. Sample SECONDARY PARTICLES from secondary_particles groups + # =========================================================================== + # TODO: Add secondary particle sampling + + +# No fission for protons + + +# ====================================================================================== +# Misc +# ====================================================================================== + +@njit +def sample_mcs_angle(E, distance, density, X0): + sigma = highland_lynch_dahl_sigma(E, distance, density, X0) + + if sigma < 0.0: + raise ValueError(f'negative sigma = {sigma}') + + # Sample theta from the Highland distribution; phi uniformly from (0, 2pi) + theta = np.abs(np.random.normal(0, sigma)) + phi = np.random.uniform(0, 2*np.pi) + + return phi, theta + + +@njit +def highland_lynch_dahl_sigma(E, distance, density, X0): + p = np.sqrt(E * (E + 2.0 * PROTON_MASS)) + beta = p / (E + PROTON_MASS) + z = 1 # Incident particle is a proton, Z=1 + + # X0 is measured in g/cm^2 + # Highland formula, modified by Lynch & Dahl + radiation_distance_fraction = density * distance / X0 + sigma = (13.6e6 / p*beta) * z * np.sqrt(radiation_distance_fraction) * (1 + 0.088 * np.log10(radiation_distance_fraction)) + sigma = np.abs(sigma) + + if sigma < 0.0: + print(f'radiation_distance_fraction = {radiation_distance_fraction}') + print(f'p = {p}, beta = {beta}, z = {z}') + print(f'density = {density}, distance = {distance}') + raise ValueError(f"negative sigma = {sigma}") + + return sigma + +@njit +def rotate_direction(particle, phi, theta): + """ + Rotate direction vector (ux, uy, uz) by polar angle theta + and azimuthal angle phi in the local frame. + Returns new (ux, uy, uz). + """ + + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + + sin_theta = np.sin(theta) + cos_theta = np.cos(theta) + cos_phi = np.cos(phi) + sin_phi = np.sin(phi) + + # Build local perpendicular axes + d = np.array([ux, uy, uz]) + perp = np.array([1.0, 0.0, 0.0]) if abs(ux) < 0.9 else np.array([0.0, 1.0, 0.0]) + u = np.cross(d, perp); u /= np.linalg.norm(u) + v = np.cross(d, u) + + d_new = (cos_theta * d + + sin_theta * cos_phi * u + + sin_theta * sin_phi * v) + d_new /= np.linalg.norm(d_new) + + particle["ux"] = d_new[0] + particle["uy"] = d_new[1] + particle["uz"] = d_new[2] + + +@njit +def calculate_total_stopping_power(particle_container, simulation, data): + particle = particle_container[0] + material = simulation["native_materials"][particle["material_ID"]] + E = particle["E"] + + total_stopping_power = 0.0 + total_rho_gcm3 = 0.0 + total_Z = 0.0 + total_A = 0.0 + # Find the total stopping power by summing over every nuclide in the material + for i in range(material["N_nuclide"]): + nuclide_ID = int(mcdc_get.native_material.nuclide_IDs(i, material, data)) + nuclide = simulation["nuclides"][nuclide_ID] + + # If no stopping power provided, we calculate it ourselves here + if not material["stopping_power_provided"]: + dedx_values = mcdc_get.nuclide.stopping_power_all(nuclide, data) + dedx_energies = mcdc_get.nuclide.stopping_power_energy_grid_all(nuclide, data) + + # TODO: replace np.interp with a non-numpy function?? + dedx = np.interp(E / 1e6, dedx_energies, dedx_values) + total_stopping_power += dedx * 1e6 + + # Convert atoms/barn-cm to g/cm3: + atomic_mass = nuclide["atomic_weight_ratio"] # mass in amu + nuclide_density = mcdc_get.native_material.nuclide_densities(i, material, data) + density_gcm3 = nuclide_density * 1e24 * atomic_mass / (6.022e23) + total_rho_gcm3 += density_gcm3 + + total_Z += nuclide["atomic_number"] + total_A += nuclide["mass_number"] + + average_Z = total_Z / material["N_nuclide"] + average_A = total_A / material["N_nuclide"] + + if material["stopping_power_provided"]: + dedx_values = mcdc_get.native_material.stopping_power_all(material, data) + dedx_energies = mcdc_get.native_material.stopping_power_energy_grid_all(material, data) + + dedx = np.interp(E / 1e6, dedx_energies, dedx_values) + total_stopping_power = dedx * 1e6 + + return average_A, average_Z, total_stopping_power, total_rho_gcm3 \ No newline at end of file diff --git a/mcdc/transport/physics/util.py b/mcdc/transport/physics/util.py index 3475a1510..8788aefce 100644 --- a/mcdc/transport/physics/util.py +++ b/mcdc/transport/physics/util.py @@ -31,6 +31,18 @@ def evaluate_electron_xs_energy_grid(e, element, data): return idx, e0, e1 +@njit +def evaluate_proton_xs_energy_grid(e, nuclide, data): + offset = nuclide["proton_xs_energy_grid_offset"] + length = nuclide["proton_xs_energy_grid_length"] + energy_grid = data[offset : offset + length] + + idx = find_bin(e, energy_grid) + e0 = energy_grid[idx] + e1 = energy_grid[idx + 1] + return idx, e0, e1 + + @njit def scatter_direction(ux, uy, uz, mu0, azi): cos_azi = math.cos(azi) diff --git a/mcdc/transport/simulation.py b/mcdc/transport/simulation.py index aca16d88f..c0c998f99 100644 --- a/mcdc/transport/simulation.py +++ b/mcdc/transport/simulation.py @@ -324,6 +324,10 @@ def step_particle(particle_container, program, data): if particle["event"] & EVENT_TIME_BOUNDARY: particle["alive"] = False + # CSDA energy depostiion + if particle["event"] & EVENT_CSDA_EDEP: + pass + # ================================================================================== # Apply techniques # ================================================================================== @@ -340,7 +344,6 @@ def step_particle(particle_container, program, data): if simulation["global_weight_roulette"]["active"]: technique.global_weight_roulette(particle_container, simulation) - @njit def move_to_event(particle_container, simulation, data): settings = simulation["settings"] @@ -393,6 +396,10 @@ def move_to_event(particle_container, simulation, data): # Distance to next collision d_collision = physics.collision_distance(particle_container, simulation, data) + # Distance to max energy loss as dictated by CSDA + if settings["csda"]: + d_csda = physics.csda_distance(particle_container, simulation, data) + # ================================================================================== # Determine event(s) # ================================================================================== @@ -422,6 +429,20 @@ def move_to_event(particle_container, simulation, data): particle["event"] = EVENT_TIME_BOUNDARY particle["surface_ID"] = -1 + # Check distance to max energy loss from CSDA + if settings["csda"]: + if d_csda < distance - COINCIDENCE_TOLERANCE: + distance = d_csda + particle["event"] = EVENT_CSDA_EDEP + particle["surface_ID"] = -1 + elif geometry.check_coincidence(d_csda, distance): + particle["event"] += EVENT_CSDA_EDEP + + if distance < 0.0: + print(f'distance = {distance}') + print(f'd_coll = {d_collision}, d_csda = {d_csda}, d_bnd = {d_boundary}') + raise ValueError(f"Negative distance") + # ================================================================================== # Move particle # ================================================================================== @@ -444,6 +465,27 @@ def move_to_event(particle_container, simulation, data): # Move particle particle_module.move(particle_container, distance, simulation, data) + # CSDA calculates energy loss after particle has moved + if settings["csda"]: + collision_data_container = np.zeros(1, type_.collision_data) + physics.csda_edep( + particle_container, collision_data_container, distance, simulation, data + ) + + # Score collision tallies (edep is a collision tally) + # TODO: maybe make edep a potential tracklength tally for CSDA? + if simulation["cycle_active"]: + cell = simulation["cells"][particle["cell_ID"]] + for i in range(cell["N_collision_tally"]): + tally_ID = int(mcdc_get.cell.collision_tally_IDs(i, cell, data)) + tally = simulation["collision_tallies"][tally_ID] + tally_module.score.collision_tally( + particle_container, + collision_data_container, + tally, + simulation, + data, + ) @njit def surface_crossing(P_arr, simulation, data): diff --git a/mcdc/transport/tally/closeout.py b/mcdc/transport/tally/closeout.py index 96f2cd66c..a6a29cbac 100644 --- a/mcdc/transport/tally/closeout.py +++ b/mcdc/transport/tally/closeout.py @@ -137,7 +137,9 @@ def _finalize(tally, simulation, data): # Check for round-off error (TODO: Check why this is needed.) if abs(radicand) < 1e-16: data[offset_sum_square + i] = 0.0 - else: + if radicand < 0.0 and abs(radicand) < 1e-6: + data[offset_sum_square + i] = 0.0 + else: data[offset_sum_square + i] = math.sqrt(radicand) diff --git a/tools/data_library_generator/neutron/generate.py b/tools/data_library_generator/neutron/generate.py index 809d12108..eacc4b9da 100644 --- a/tools/data_library_generator/neutron/generate.py +++ b/tools/data_library_generator/neutron/generate.py @@ -1,8 +1,8 @@ -import ACEtk import argparse import h5py import numpy as np import os +import ACEtk from tqdm import tqdm diff --git a/tools/data_library_generator/neutron/util.py b/tools/data_library_generator/neutron/util.py index 44f9b5897..9a83892e3 100644 --- a/tools/data_library_generator/neutron/util.py +++ b/tools/data_library_generator/neutron/util.py @@ -45,7 +45,16 @@ def decode_ace_name(name: str): S = offset // 100 A = offset % 100 - T = ACE_TEMPERATURE_LIB81[extension] + # Proton data: ENDF70PROT + if extension == "70h": + T = 293.6 + + # Proton data: TENDL-19 (defaults at 0K, I think) + if extension == "19h": + T = 0 + + else: + T = ACE_TEMPERATURE_LIB81[extension] return Z, A, S, T diff --git a/tools/data_library_generator/proton/README.md b/tools/data_library_generator/proton/README.md new file mode 100644 index 000000000..fd4847f8e --- /dev/null +++ b/tools/data_library_generator/proton/README.md @@ -0,0 +1,167 @@ +# MC/DC Proton Data Library Generator +Converts ACE-format proton cross section data from TENDL2021 into MC/DC's +per-nuclide HDF5 format for proton transport. Also uses the NIST PSTAR +database to add stopping power data to the HDF5 files. + +## Prerequisites + +- Installing ACEtk from source: [link](https://github.com/njoy/ACEtk) +- Dependencies: `pip install h5py numpy tqdm` +- For any nuclide that you want stopping power for, download the stopping power data from NIST'S PSTAR database: https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html + - When downloading the stopping power data, select only the total stopping power. This generator script uses the load_pstar_file function to extract the data assuming two columns of data. +- You need the TENDL2021 ACE files for protons. Avaiable at this link: https://tendl.imperial.ac.uk/tendl_2021/tar.html + - Or, download the tar file directly: https://tendl.imperial.ac.uk/tendl_2021/tar_files/TENDL-ACE-p.tgz + +## Environment Variables +| Variable | Description | +|------------------------|------------------------------------------------------------------------| +| `MCDC_ACELIB_PROTON` | Path to the directory containing the TENDL2021 ACE files. | +| `MCDC_LIB_PROTON` | Path to the output directory for MC/DC HDF5 files. | +| `MCDC_PSTAR_LIB` | Path to the directory containing the PSTAR stopping power table files. | + +## Usage +```bash +export MCDC_ACELIB_PROTON=/path/to/tendl2021/acefiles +export MCDC_LIB_PROTON=/path/to/mcdc/proton/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=103) in the TENDL2021 ACE file library, the generator: +1. Loads the data from the ACE table, and writes basic data (name, temperature, mass, etc.) to the HDF5 file. +2. Extracts the stopping power data (if present in $MCDC_PSTAR_LIB). +3. Extracts the MT numbers for elastic scattering reactions and their cross sections. +4. Extracts the MT numbers for capture reactions and their cross sections. +5. Extracts the MT numbers for inelastic scattering reactions and their cross sections. +6. Extracts the energy & angular distributions for scattering reactions. +7. Creates a data block in the HDF5 file to handle the secondary particle products & energies. +8. If there are isotopes present with PSTAR stopping power data, but without an ACE file from TENDL, it creates + an HDF5 file that contains the stopping power data. + +## Output HDF5 Schema +``` +File attrs: source_title, source_version, source_date, source_comments (if present) + +-K.h5 +├── nuclide_name (string) +├── excitation_level (int) +├── temperature (float; attr: unit="K") +├── atomic_number (int) +├── mass_number (int) +├── atomic_weight_ratio (float) +├── radiation_length (float; attr: unit="g/cm2") +├── fissionable (bool) +├── stopping_power/ (present only if a matching PSTAR file was found) +│ ├── energy (1-D array; attr: unit="MeV") +│ └── total_stopping_power (1-D array; attr: unit="MeV cm2/g") +├── proton_reactions/ +│ ├── xs_energy_grid (1-D array; attr: unit="MeV") +│ ├── elastic_scattering/ +│ │ └── MT-002/ (attr: MT=2) +│ │ ├── xs (1-D array, barns; attr: offset=0) +│ │ ├── Q-value (float=0.0; attr: unit="MeV") +│ │ ├── reference_frame (string: "COM") +│ │ └── angular_cosine_distribution/ (attr: type="energy-correlated"; see [A] below) +│ ├── capture/ (one group per capture MT, i.e. multiplicity=0) +│ │ └── MT-NNN/ (attr: MT) +│ │ ├── xs (1-D array, barns; attr: offset) +│ │ ├── Q-value (float; attr: unit="MeV") +│ │ └── reference_frame (string: "LAB" or "COM") +│ ├── inelastic_reaction/ (present only if any inelastic MTs exist) +│ │ └── MT-NNN/ (attr: MT; one per inelastic MT) +│ │ ├── xs (1-D array, barns; attr: offset) +│ │ ├── Q-value (float; attr: unit="MeV") +│ │ ├── reference_frame (string: "LAB" or "COM") +│ │ ├── multiplicity (int) +│ │ ├── angular_cosine_distribution/ (see [A] below) +│ │ ├── spectrum_probability_grid (1-D array; attr: unit="MeV") +│ │ ├── spectrum_probability (2-D array [grid x n_dist]) +│ │ └── energy_spectrum-N/ (one per outgoing-energy law; see [B] below) +│ └── fission/ (present only if fissionable) +│ ├── MT-NNN/ (attr: MT; MT-018 or the fission-chance MTs 19/20/21/38) +│ │ ├── xs (1-D array, barns; attr: offset) +│ │ ├── Q-value (float; attr: unit="MeV") +│ │ ├── reference_frame (string: "LAB" or "COM") +│ │ ├── angular_cosine_distribution/ (see [A] below) +│ │ ├── spectrum_probability_grid (1-D array; attr: unit="MeV") +│ │ ├── spectrum_probability (2-D array [grid x n_dist]) +│ │ └── energy_spectrum-N/ (see [B] below) +│ ├── prompt_multiplicity/ (see [C] below) +│ ├── delayed_multiplicity/ (optional; see [C] below) +│ └── delayed_neutron_precursors/ (optional) +│ ├── fractions (1-D array, one per precursor group) +│ ├── decay_rates (1-D array; attr: unit="/s") +│ └── energy_spectrum-N/ (one per precursor group; see [B] below) +└── secondary_particles/ (present only if ACE table has secondary-particle data) + └── ZAP_/ (attrs: ZAP, particle_name; one per secondary particle type) + └── MT-NNN/ (attrs: MT, multiplicity, reference_frame) + ├── production_xs (1-D array, barns; attr: offset) + ├── kalbach_mann/ (attr: type="kalbach-mann") + │ ├── energy (1-D array; attr: unit="MeV") + │ ├── offset (1-D int array, one entry per incident energy) + │ ├── energy_out (1-D array; attr: unit="MeV") + │ ├── pdf (1-D array) + │ ├── cdf (1-D array) + │ ├── precompound_factor (1-D array) + │ └── angular_slope (1-D array) + └── angular_cosine_distribution/ (see [A] below) + + +[A] angular_cosine_distribution/ schema (load_cosine_distribution): + attr: type="given_in_energy_distribution" (angular data embedded in the energy-distribution block — no other content) + — or — + attr: type="tabulated"; attr: unit="MeV" + ├── incident_energies (1-D array) + └── E_in_i/ (one group per incident energy) + attr: type="tabulated" → cosines, pdf, cdf (1-D arrays) + attr: type="isotropic" → (no datasets) + +[B] energy_spectrum-N/ schema (load_energy_distribution), attr: law = ENDF law number: + law=44 (Kalbach-Mann): + attr: type="kalbach-mann" + ├── energy (1-D array; attr: unit="MeV") + ├── offset (1-D int array) + ├── energy_out (1-D array; attr: unit="MeV") + ├── pdf, cdf (1-D arrays) + ├── precompound_factor (1-D array) + └── angular_slope (1-D array) + law=4 (tabulated outgoing energy): + ├── incident_energies (1-D array) + └── E_in_k/ → outgoing_energies, pdf, cdf (1-D arrays) + law=3 (level scattering): + ├── C1 (float) + └── C2 (float) + law=1 (equiprobable bins): + ├── incident_energies (1-D array) + └── E_in_k/ → energies (1-D array) + law=-1 (unrecognized type): + attr: type_name= + └── xss_array (1-D array; only if extraction succeeds) + +[C] prompt_multiplicity/ and delayed_multiplicity/ schema (load_fission_multiplicity): + attr: type="tabulated" → energies, multiplicities (1-D arrays) + attr: type="polynomial" → coefficients (1-D array) + attr: type="unknown" → attr: type_name= (no data) + + +── Stopping-power-only fallback (process_pstar_only_file; H-1, H-2, He-3, He-4 when no ACE file exists) ── + +-K.h5 +├── nuclide_name (string) +├── excitation_level (int = 0) +├── temperature (float; attr: unit="K") +├── atomic_number (int) +├── mass_number (int) +├── atomic_weight_ratio (float) +├── radiation_length (float; attr: unit="g/cm2") +├── fissionable (bool = False) +└── stopping_power/ + ├── energy (1-D array; attr: unit="MeV") + └── total_stopping_power (1-D array; attr: unit="MeV cm2/g") +``` + +## See Also +- [TENDL2021](https://tendl.imperial.ac.uk/tendl_2021/tendl2021.html) diff --git a/tools/data_library_generator/proton/generate.py b/tools/data_library_generator/proton/generate.py new file mode 100644 index 000000000..d51fa456c --- /dev/null +++ b/tools/data_library_generator/proton/generate.py @@ -0,0 +1,822 @@ +# The majority of this script was written by Anthropic's Claude + +import argparse +import os +import sys + +import h5py +import numpy as np +from tqdm import tqdm +import ACEtk + + +# -- Constants ----------------------------------------------------------------- + +ZAP_NAMES = { + 1: "neutron", + 31: "deuteron", + 32: "triton", + 33: "He3", + 34: "alpha", +} + +Z_TO_SYMBOL = { + 1:"H", 2:"He", 3:"Li", 4:"Be", 5:"B", 6:"C", 7:"N", 8:"O", + 9:"F", 10:"Ne", 11:"Na", 12:"Mg", 13:"Al", 14:"Si", 15:"P", 16:"S", + 17:"Cl", 18:"Ar", 19:"K", 20:"Ca", 21:"Sc", 22:"Ti", 23:"V", 24:"Cr", + 25:"Mn", 26:"Fe", 27:"Co", 28:"Ni", 29:"Cu", 30:"Zn", 31:"Ga", 32:"Ge", + 33:"As", 34:"Se", 35:"Br", 36:"Kr", 37:"Rb", 38:"Sr", 39:"Y", 40:"Zr", + 41:"Nb", 42:"Mo", 43:"Tc", 44:"Ru", 45:"Rh", 46:"Pd", 47:"Ag", 48:"Cd", + 49:"In", 50:"Sn", 51:"Sb", 52:"Te", 53:"I", 54:"Xe", 55:"Cs", 56:"Ba", + 57:"La", 58:"Ce", 59:"Pr", 60:"Nd", 61:"Pm", 62:"Sm", 63:"Eu", 64:"Gd", + 65:"Tb", 66:"Dy", 67:"Ho", 68:"Er", 69:"Tm", 70:"Yb", 71:"Lu", 72:"Hf", + 73:"Ta", 74:"W", 75:"Re", 76:"Os", 77:"Ir", 78:"Pt", 79:"Au", 80:"Hg", + 81:"Tl", 82:"Pb", 83:"Bi", 84:"Po", 85:"At", 86:"Rn", 87:"Fr", 88:"Ra", + 89:"Ac", 90:"Th", 91:"Pa", 92:"U", 93:"Np", 94:"Pu", 95:"Am", 96:"Cm", + 97:"Bk", 98:"Cf", 99:"Es",100:"Fm",101:"Md",102:"No",103:"Lr",104:"Rf", + 105:"Db",106:"Sg",107:"Bh",108:"Hs",109:"Mt",110:"Ds",111:"Rg",112:"Cn", + 113:"Nh",114:"Fl",115:"Mc",116:"Lv",117:"Ts",118:"Og" +} + +RADIATION_LENGTH_FROM_Z = { + "H": 63.04, "He": 94.32, "Li": 82.77, "Be": 65.19, "B": 52.68, "C": 42.70, + "N": 37.99, "O": 34.24, "F": 32.93, "Ne": 28.93, "Na": 27.74, "Mg": 25.03, + "Al": 24.01, "Si": 21.82, "P": 21.21, "S": 19.50, "Cl": 19.28, "Ar": 19.55, + "K": 17.32, "Ca": 16.14, "Sc": 16.55, "Ti": 16.16, "V": 15.84, "Cr": 14.94, + "Mn": 14.64, "Fe": 13.84, "Co": 13.62, "Ni": 12.68, "Cu": 12.86, "Zn": 12.43, + "Ga": 12.47, "Ge": 12.25, "As": 11.94, "Se": 11.91, "Br": 11.42, "Kr": 11.37, + "Rb": 11.03, "Sr": 10.76, "Y": 10.41, "Zr": 10.20, "Nb": 9.92, "Mo": 9.80, + "Tc": 9.58, "Ru": 9.48, "Rh": 9.27, "Pd": 9.20, "Ag": 8.97, "Cd": 9.00, + "In": 8.85, "Sn": 8.82, "Sb": 8.73, "Te": 8.83, "I": 8.48, "Xe": 8.48, + "Cs": 8.31, "Ba": 8.31, "La": 8.14, "Ce": 7.96, "Pr": 7.76, "Nd": 7.71, + "Pm": 7.51, "Sm": 7.57, "Eu": 7.44, "Gd": 7.48, "Tb": 7.36, "Dy": 7.32, + "Ho": 7.23, "Er": 7.14, "Tm": 7.03, "Yb": 7.02, "Lu": 6.92, "Hf": 6.89, + "Ta": 6.82, "W": 6.76, "Re": 6.69, "Os": 6.68, "Ir": 6.59, "Pt": 6.54, + "Au": 6.46, "Hg": 6.44, "Tl": 6.42, "Pb": 6.37, "Bi": 6.29, "Po": 6.16, + "At": 6.07, "Rn": 6.28, "Fr": 6.19, "Ra": 6.15, "Ac": 6.06, "Th": 6.07, + "Pa": 5.93, "U": 6.00, "Np": 5.87, "Pu": 5.93, "Am": 5.80, "Cm": 5.79, + "Bk": 5.69, "Cf": 5.68, "Es": 5.61, "Fm": 5.62, "Md": 5.55, "No": 5.48, + "Lr": 5.45, "Rf": 5.47, "Db": 5.40, "Sg": 5.34, "Bh": 5.27, "Hs": 5.17, + "Mt": 5.26, "Ds": 5.24, "Rg": 5.18, "Cn": 5.16, "Nh": 5.10, "Fl": 5.08, + "Mc": 5.01, "Lv": 5.00, "Ts": 4.95, "Og": 4.88, +} + +SYMBOL_TO_Z = {v: k for k, v in Z_TO_SYMBOL.items()} + +# Isotopes to generate stopping-power-only HDF5 files for when no ACE file +# exists. Covers H and He which TENDL excludes because TALYS doesn't apply. +# Format: (symbol, A, atomic_weight_ratio) +# AWR = atomic mass / neutron mass; neutron mass = 1.008664916 u +PSTAR_ONLY_ISOTOPES = [ + ("H", 1, 1.00794 / 1.008664916), # natural H ≈ H-1 + ("H", 2, 2.01410 / 1.008664916), # deuterium + ("He", 3, 3.01603 / 1.008664916), # He-3 + ("He", 4, 4.00260 / 1.008664916), # He-4 +] + +# Redundant sum MTs that must not be double-counted +REDUNDANT_MTS = [1, 3, 4, 10, 101, 103, 104, 105, 106, 107] +FISSION_CHANCE_MTS = [19, 20, 21, 38] + +# Temperature written to all files (TENDL proton ACE files report 0 K) +T_KELVIN = 0.0 + + +# -- Utility ------------------------------------------------------------------- + +def print_error(msg): + print(f"\n[ERROR] {msg}", file=sys.stderr) + raise ValueError(msg) + sys.exit(1) + + +def print_note(msg): + print(f" [note] {msg}") + + +def decode_ace_zaid(zaid): + """Return (Z, A, S, T=0) from an ACE ZAID string.""" + za = int(zaid.strip().split(".")[0]) + S = 0 + if za >= 600000: + S = (za % 1000) // 400 + za = za - S * 400 + return za // 1000, za % 1000, S, 0 + + +def load_pstar_file(filepath): + """ + Load a two-column PSTAR file (Energy MeV, Stopping power MeV cm^2/g). + Returns (energies, stopping_powers) as float64 arrays. + """ + energies, sps = [], [] + with open(filepath) as f: + for line in f: + parts = line.strip().split() + if len(parts) != 2: + continue + try: + energies.append(float(parts[0])) + sps.append(float(parts[1])) + except ValueError: + continue + return np.array(energies), np.array(sps) + + +def write_stopping_power(file, pstar_dir, symbol, verbose=False): + """ + Write stopping_power group into an open HDF5 file if a PSTAR file exists. + Returns True if data was written. + """ + if pstar_dir is None: + return False + pstar_path = os.path.join(pstar_dir, f"{symbol}.txt") + if not os.path.exists(pstar_path): + if verbose: + print(f" [warn] No PSTAR file for {symbol}") + return False + if verbose: + print(f" Loading PSTAR from {pstar_path}") + E_s, S_s = load_pstar_file(pstar_path) + sp = file.create_group("stopping_power") + sp.create_dataset("energy", data=E_s).attrs["unit"] = "MeV" + sp.create_dataset("total_stopping_power", data=S_s).attrs["unit"] = "MeV cm2/g" + return True + + +# -- Distribution writers ------------------------------------------------------ + +def load_cosine_distribution(data, h5_group): + """ + Write a tabulated angular distribution into h5_group. + Returns False if the distribution is embedded in a Kalbach-Mann block + (DistributionGivenElsewhere), True otherwise. + """ + if isinstance(data, ACEtk.continuous.DistributionGivenElsewhere): + h5_group.attrs["type"] = "given_in_energy_distribution" + return False + + h5_group.attrs["type"] = "tabulated" + h5_group.attrs["unit"] = "MeV" + h5_group.create_dataset("incident_energies", data=np.array(data.incident_energies)) + + for i, subdist in enumerate(data.distributions): + eg = h5_group.create_group(f"E_in_{i + 1}") + if isinstance(subdist, ACEtk.continuous.TabulatedAngularDistribution): + eg.attrs["type"] = "tabulated" + eg.create_dataset("cosines", data=np.array(subdist.cosines)) + eg.create_dataset("pdf", data=np.array(subdist.pdf)) + eg.create_dataset("cdf", data=np.array(subdist.cdf)) + else: + eg.attrs["type"] = "isotropic" + + return True + + +def _write_kalbach_mann(km_data, h5_group): + """ + Write a KalbachMannDistributionData into h5_group as flat arrays. + offset[i] gives the starting index in the flat arrays for incident energy i. + """ + h5_group.attrs["type"] = "kalbach-mann" + + NE = km_data.number_incident_energies + h5_group.create_dataset( + "energy", data=np.array(km_data.incident_energies) + ).attrs["unit"] = "MeV" + + offset, energy_out, pdf, cdf, r_vals, a_vals = [], [], [], [], [], [] + for i in range(1, NE + 1): + dist = km_data.distribution(i) + offset.append(len(energy_out)) + energy_out.extend(dist.outgoing_energies) + pdf.extend(dist.pdf) + cdf.extend(dist.cdf) + r_vals.extend(dist.precompound_fraction_values) + a_vals.extend(dist.angular_distribution_slope_values) + + h5_group.create_dataset("offset", data=np.array(offset, dtype=np.int32)) + h5_group.create_dataset( + "energy_out", data=np.array(energy_out) + ).attrs["unit"] = "MeV" + h5_group.create_dataset("pdf", data=np.array(pdf)) + h5_group.create_dataset("cdf", data=np.array(cdf)) + h5_group.create_dataset("precompound_factor", data=np.array(r_vals)) + h5_group.create_dataset("angular_slope", data=np.array(a_vals)) + + +def load_energy_distribution(data, h5_group): + """Write a primary-particle outgoing energy distribution into h5_group.""" + if isinstance(data, ACEtk.continuous.KalbachMannDistributionData): + h5_group.attrs["law"] = 44 + _write_kalbach_mann(data, h5_group) + + elif isinstance(data, ACEtk.continuous.OutgoingEnergyDistributionData): + h5_group.attrs["law"] = 4 + h5_group.create_dataset( + "incident_energies", data=np.array(data.incident_energies) + ) + for k, dist in enumerate(data.distributions): + eg = h5_group.create_group(f"E_in_{k + 1}") + eg.create_dataset("outgoing_energies", data=np.array(dist.outgoing_energies)) + eg.create_dataset("pdf", data=np.array(dist.pdf)) + eg.create_dataset("cdf", data=np.array(dist.cdf)) + + elif isinstance(data, ACEtk.continuous.LevelScatteringData): + h5_group.attrs["law"] = 3 + h5_group.create_dataset("C1", data=data.C1) + h5_group.create_dataset("C2", data=data.C2) + + elif isinstance(data, ACEtk.continuous.EquiprobableOutgoingEnergyBins): + h5_group.attrs["law"] = 1 + h5_group.create_dataset( + "incident_energies", data=np.array(data.incident_energies) + ) + for k, dist in enumerate(data.distributions): + h5_group.create_group(f"E_in_{k + 1}").create_dataset( + "energies", data=np.array(dist.energies) + ) + + else: + h5_group.attrs["law"] = -1 + h5_group.attrs["type_name"] = type(data).__name__ + try: + h5_group.create_dataset("xss_array", data=np.array(data.xss_array)) + except Exception: + pass + + +def load_fission_multiplicity(data, h5_group): + if isinstance(data, ACEtk.continuous.TabulatedFissionMultiplicity): + h5_group.attrs["type"] = "tabulated" + h5_group.create_dataset("energies", data=np.array(data.energies)) + h5_group.create_dataset("multiplicities", data=np.array(data.multiplicities)) + elif isinstance(data, ACEtk.continuous.PolynomialFissionMultiplicity): + h5_group.attrs["type"] = "polynomial" + h5_group.create_dataset("coefficients", data=np.array(data.coefficients)) + else: + h5_group.attrs["type"] = "unknown" + h5_group.attrs["type_name"] = type(data).__name__ + + +# -- Secondary particles ------------------------------------------------------- + +def load_secondary_particles(ace_table, file, verbose=False): + n_types = ace_table.number_secondary_particle_types + if n_types == 0: + return + + type_block = ace_table.secondary_particle_type_block + info_block = ace_table.secondary_particle_information_block + rx_block = ace_table.secondary_particle_reaction_number_block + tyr_block = ace_table.secondary_particle_frame_and_multiplicity_block + xs_block = ace_table.secondary_particle_production_cross_section_block + edy_block = ace_table.secondary_particle_energy_distribution_block + + has_ang = False + try: + ang_block = ace_table.secondary_particle_angular_distribution_block + has_ang = True + except Exception: + pass + + sec_group = file.create_group("secondary_particles") + + pi_method = next( + (c for c in ["particle_identifier", "ZAP", "type", "particle_type"] + if hasattr(type_block, c)), + None + ) + if pi_method is None: + raise AttributeError( + f"Cannot find particle identifier on {type(type_block).__name__}. " + f"Available: {[x for x in dir(type_block) if not x.startswith('_')]}" + ) + + for i in range(1, n_types + 1): + zap = getattr(type_block, pi_method)(i) + name = ZAP_NAMES.get(zap, f"ZAP_{zap}") + n_rx = int(info_block.number_reactions[i - 1]) + + if verbose: + print(f" Secondary type {i}: ZAP={zap} ({name}), {n_rx} reactions") + + zap_group = sec_group.create_group(f"ZAP_{zap}") + zap_group.attrs["ZAP"] = zap + zap_group.attrs["particle_name"] = name + + rx_i = rx_block(i) + tyr_i = tyr_block(i) + xs_i = xs_block(i) + edy_i = edy_block(i) + ang_i = ang_block(i) if has_ang else None + + xs_method = next( + (c for c in ["cross_sections", "cross_section", "xs"] if hasattr(xs_i, c)), + None + ) + off_method = next( + (c for c in ["energy_index", "offset", "locator", "index"] if hasattr(xs_i, c)), + None + ) + edy_method = next( + (c for c in ["energy_distribution_data", "distribution_data", "distribution"] + if hasattr(edy_i, c)), + None + ) + + for j in range(1, n_rx + 1): + MT = rx_i.MT(j) + nu_raw = tyr_i.multiplicity(j) + nu = nu_raw - 100 if nu_raw >= 100 else nu_raw + rf_raw = tyr_i.reference_frame(j) + rf = ("LAB" if rf_raw == ACEtk.ReferenceFrame.Laboratory else + "COM" if rf_raw == ACEtk.ReferenceFrame.CentreOfMass else str(rf_raw)) + + mt = zap_group.create_group(f"MT-{MT:03}") + mt.attrs["MT"] = MT + mt.attrs["multiplicity"] = nu + mt.attrs["reference_frame"] = rf + + if verbose: + print(f" MT={MT:03} nu_raw={nu_raw} nu={nu} frame={rf}") + + empty_xs = np.zeros(0, dtype=float) + if xs_method and off_method: + try: + ds = mt.create_dataset( + "production_xs", data=np.array(getattr(xs_i, xs_method)(j)) + ) + ds.attrs["offset"] = int(getattr(xs_i, off_method)(j)) - 1 + ds.attrs["unit"] = "barns" + except Exception as exc: + ds = mt.create_dataset("production_xs", data=empty_xs) + ds.attrs["offset"] = 0 + ds.attrs["unit"] = "barns" + if verbose: + print(f" [warn] production xs: {exc}") + else: + ds = mt.create_dataset("production_xs", data=empty_xs) + ds.attrs["offset"] = 0 + ds.attrs["unit"] = "barns" + if verbose: + print(f" [warn] xs methods not found: " + f"{[x for x in dir(xs_i) if not x.startswith('_')]}") + + if edy_method: + try: + _write_kalbach_mann( + getattr(edy_i, edy_method)(j), + mt.create_group("kalbach_mann") + ) + except Exception as exc: + if verbose: + print(f" [warn] energy dist: {exc}") + elif verbose: + print(f" [warn] edy method not found: " + f"{[x for x in dir(edy_i) if not x.startswith('_')]}") + + if ang_i is not None: + try: + load_cosine_distribution( + ang_i.angular_distribution_data(j), + mt.create_group("angular_cosine_distribution") + ) + except Exception: + pass + + +# -- Per-file processing ------------------------------------------------------- + +def process_ace_file(ace_path, output_dir, pstar_dir=None, verbose=False): + """Convert a single ACE proton file to HDF5. Returns the output filename.""" + with open(ace_path) as f: + header = ACEtk.Header.from_string(f.readline()) + + Z, A, S, _ = decode_ace_zaid(header.zaid) + symbol = Z_TO_SYMBOL.get(Z, f"Z{Z}") + nuclide_name = f"{symbol}{A}" if S == 0 else f"{symbol}{A}m{S}" + + # Special case for deuterium + if symbol == "H2": + radiation_length = 125.98 + else: + radiation_length = RADIATION_LENGTH_FROM_Z.get(symbol) + + ace_table = ACEtk.ContinuousEnergyTable.from_file(ace_path) + mcdc_name = f"{nuclide_name}-{T_KELVIN}K.h5" + out_path = os.path.join(output_dir, mcdc_name) + + if verbose: + print(f"\n{'='*80}") + print(f" {os.path.basename(ace_path)} -> {mcdc_name}") + print(f" Z={Z} A={A} S={S} T={T_KELVIN} K") + + file = h5py.File(out_path, "w") + + # Metadata + hdr = ace_table.header + file.attrs["source_title"] = hdr.title + file.attrs["source_version"] = hdr.version + file.attrs["source_date"] = hdr.date + if hasattr(hdr, "comments"): + file.attrs["source_comments"] = hdr.comments + + file.create_dataset("nuclide_name", data=nuclide_name) + file.create_dataset("excitation_level", data=S) + file.create_dataset("temperature", data=T_KELVIN).attrs["unit"] = "K" + file.create_dataset("atomic_number", data=ace_table.atom_number) + file.create_dataset("mass_number", data=ace_table.mass_number) + file.create_dataset("atomic_weight_ratio", data=ace_table.atomic_weight_ratio) + file.create_dataset("radiation_length", data=radiation_length).attrs["unit"] = "g/cm2" + fissionable = ace_table.fission_multiplicity_block is not None + file.create_dataset("fissionable", data=fissionable) + + write_stopping_power(file, pstar_dir, symbol, verbose=verbose) + + # Reaction classification + nu_block = ace_table.frame_and_multiplicity_block + rx_block = ace_table.reaction_number_block + N_reaction = nu_block.number_reactions + + proton_reactions = file.create_group("proton_reactions") + elastic_group = proton_reactions.create_group("elastic_scattering") + capture_group = proton_reactions.create_group("capture") + inelastic_group = proton_reactions.create_group("inelastic_scattering") + fission_group = proton_reactions.create_group("fission") + + elastic_MTs = [2] + capture_MTs = [] + inelastic_MTs = [] + fission_MTs = ([18] if rx_block.has_MT(18) else + [MT for MT in FISSION_CHANCE_MTS if rx_block.has_MT(MT)]) + + for i in range(N_reaction): + idx = i + 1 + MT = rx_block.MT(idx) + if MT in REDUNDANT_MTS + elastic_MTs + fission_MTs or MT > 891: + continue + nu_raw = nu_block.multiplicity(idx) + if not isinstance(nu_raw, int): + print_error(f"Non-integer multiplicity for MT-{MT:03} in {ace_path}") + nu = nu_raw - 100 if nu_raw >= 100 else nu_raw + if nu == 0: capture_MTs.append(MT) + elif nu > 0: inelastic_MTs.append(MT) + else: print_error(f"Negative multiplicity for MT-{MT:03} in {ace_path}") + + for grp, mts in [(elastic_group, elastic_MTs), + (capture_group, capture_MTs), + (inelastic_group, inelastic_MTs), + (fission_group, fission_MTs)]: + for MT in mts: + grp.create_group(f"MT-{MT:03}").attrs["MT"] = MT + + if verbose: + print(f" Elastic: {elastic_MTs} Capture: {capture_MTs} " + f"Inelastic: {inelastic_MTs}" + + (f" Fission: {fission_MTs}" if fissionable else "")) + + if not fissionable: + del file["proton_reactions/fission"] + if not inelastic_MTs: + del file["proton_reactions/inelastic_scattering"] + + # Cross sections + xs0 = ace_table.principal_cross_section_block + xs_main = ace_table.cross_section_block + + proton_reactions.create_dataset( + "xs_energy_grid", data=np.array(xs0.energies) + ).attrs["unit"] = "MeV" + + ds = elastic_group.create_dataset("MT-002/xs", data=np.array(xs0.elastic)) + ds.attrs["offset"] = 0 + ds.attrs["unit"] = "barns" + + for mts, grp in [(capture_MTs, capture_group), + (inelastic_MTs, inelastic_group), + (fission_MTs, fission_group if fissionable else None)]: + if grp is None: + continue + for MT in mts: + idx = rx_block.index(MT) + ds = grp.create_dataset( + f"MT-{MT:03}/xs", data=np.array(xs_main.cross_sections(idx)) + ) + ds.attrs["offset"] = xs_main.energy_index(idx) - 1 + ds.attrs["unit"] = "barns" + + # Q-values + q_block = ace_table.reaction_qvalue_block + elastic_group.create_dataset("MT-002/Q-value", data=0.0).attrs["unit"] = "MeV" + + for mts, grp in [(capture_MTs, capture_group), + (inelastic_MTs, inelastic_group), + (fission_MTs, fission_group if fissionable else None)]: + if grp is None: + continue + for MT in mts: + idx = rx_block.index(MT) + grp.create_dataset( + f"MT-{MT:03}/Q-value", data=q_block.q_value(idx) + ).attrs["unit"] = "MeV" + + # Reference frames + elastic_group.create_dataset("MT-002/reference_frame", data="COM") + + for mts, grp in [(capture_MTs, capture_group), + (inelastic_MTs, inelastic_group), + (fission_MTs, fission_group if fissionable else None)]: + if grp is None: + continue + for MT in mts: + idx = rx_block.index(MT) + rf = nu_block.reference_frame(idx) + rf_str = ("LAB" if rf == ACEtk.ReferenceFrame.Laboratory else + "COM" if rf == ACEtk.ReferenceFrame.CentreOfMass else str(rf)) + grp.create_dataset(f"MT-{MT:03}/reference_frame", data=rf_str) + + # Inelastic multiplicities + for MT in inelastic_MTs: + idx = rx_block.index(MT) + nu_raw = nu_block.multiplicity(idx) + inelastic_group.create_dataset( + f"MT-{MT:03}/multiplicity", + data=nu_raw - 100 if nu_raw >= 100 else nu_raw + ) + + # Angular distributions + angle_block = ace_table.angular_distribution_block + + ag = elastic_group.create_group("MT-002/angular_cosine_distribution") + ag.attrs["type"] = "energy-correlated" + if not load_cosine_distribution(angle_block.angular_distribution_data(0), ag) \ + and verbose: + print_note("MT-002 angular distribution is given in energy block") + + for mts, grp in [(inelastic_MTs, inelastic_group), + (fission_MTs, fission_group if fissionable else None)]: + if grp is None: + continue + for MT in mts: + idx = rx_block.index(MT) + ag = grp.create_group(f"MT-{MT:03}/angular_cosine_distribution") + if not load_cosine_distribution( + angle_block.angular_distribution_data(idx), ag) and verbose: + print_note(f"MT-{MT:03} angular distribution is given in energy block") + + # Primary energy distributions + energy_block = ace_table.energy_distribution_block + + for mts, grp in [(inelastic_MTs, inelastic_group), + (fission_MTs, fission_group if fissionable else None)]: + if grp is None: + continue + for MT in mts: + idx = rx_block.index(MT) + data = energy_block.energy_distribution_data(idx) + + if not isinstance(data, ACEtk.continuous.MultiDistributionData): + grp.create_dataset( + f"MT-{MT:03}/spectrum_probability_grid", + data=np.array([0.0, 30.0]) + ).attrs["unit"] = "MeV" + grp.create_dataset( + f"MT-{MT:03}/spectrum_probability", data=np.array([[1.0]]) + ) + load_energy_distribution( + data, grp.create_group(f"MT-{MT:03}/energy_spectrum-1") + ) + else: + N_dist = data.number_distributions + probs = data.probabilities + + if all(p.number_interpolation_regions == 0 for p in probs): + prob_grid = np.array([0.0, 30.0]) + prob = np.zeros((1, N_dist)) + for k in range(N_dist): + prob[0, k] = max(data.probability(k + 1).probabilities) + elif (all(p.number_interpolation_regions == 1 for p in probs) + and all(p.interpolants[0] == 1 for p in probs)): + prob_grid = np.array(data.probability(1).energies) + prob = np.zeros((len(prob_grid) - 1, N_dist)) + for k in range(N_dist): + prob[:, k] = np.array(data.probability(k + 1).probabilities[:-1]) + else: + print_error( + f"Unsupported multi-distribution probability for MT-{MT:03}" + ) + + grp.create_dataset( + f"MT-{MT:03}/spectrum_probability_grid", data=prob_grid + ).attrs["unit"] = "MeV" + grp.create_dataset(f"MT-{MT:03}/spectrum_probability", data=prob) + for k in range(N_dist): + load_energy_distribution( + data.distribution(k + 1), + grp.create_group(f"MT-{MT:03}/energy_spectrum-{k + 1}") + ) + + load_secondary_particles(ace_table, file, verbose=verbose) + + # Fission data + if fissionable: + prompt_block = ace_table.fission_multiplicity_block + delayed_block = ace_table.delayed_fission_multiplicity_block + dnp_block = ace_table.delayed_neutron_precursor_block + + load_fission_multiplicity( + prompt_block.multiplicity, + fission_group.create_group("prompt_multiplicity") + ) + if delayed_block is not None: + load_fission_multiplicity( + delayed_block.multiplicity, + fission_group.create_group("delayed_multiplicity") + ) + + if dnp_block is not None: + N_DNP = dnp_block.number_delayed_precursors + fractions = np.zeros(N_DNP) + decay_rates = np.zeros(N_DNP) + for k in range(N_DNP): + d = dnp_block.precursor_group_data(k + 1) + if (d.number_interpolation_regions != 0 + or len(d.probabilities[:]) != 2 + or d.probabilities[0] != d.probabilities[1]): + print_error("Non-constant delayed neutron precursor fraction") + fractions[k] = d.probabilities[0] + decay_rates[k] = d.decay_constant + + prec = fission_group.create_group("delayed_neutron_precursors") + prec.create_dataset("fractions", data=fractions) + prec.create_dataset("decay_rates", data=decay_rates).attrs["unit"] = "/s" + + delayed_spectrum_block = ace_table.delayed_neutron_energy_distribution_block + for k in range(N_DNP): + load_energy_distribution( + delayed_spectrum_block.energy_distribution_data(k + 1), + prec.create_group(f"energy_spectrum-{k + 1}") + ) + + file.close() + return mcdc_name + + +def process_pstar_only_file(symbol, A, awr, output_dir, pstar_dir, verbose=False): + """ + Create a minimal HDF5 file for an isotope that has no ACE data but does + have a PSTAR stopping power file. Returns the output filename, or None if + no PSTAR file was found. + """ + Z = SYMBOL_TO_Z[symbol] + nuclide_name = f"{symbol}{A}" + mcdc_name = f"{nuclide_name}-{T_KELVIN}K.h5" + out_path = os.path.join(output_dir, mcdc_name) + + # Special case for deuterium + if symbol == "H2": + radiation_length = 125.98 + else: + radiation_length = RADIATION_LENGTH_FROM_Z.get(symbol) + + if verbose: + print(f"\n{'='*80}") + print(f" (no ACE) -> {mcdc_name} [stopping power only]") + + file = h5py.File(out_path, "w") + + file.attrs["source_title"] = "PSTAR (NIST) stopping power only — no ACE data" + file.attrs["source_version"] = "N/A" + file.attrs["source_date"] = "N/A" + + file.create_dataset("nuclide_name", data=nuclide_name) + file.create_dataset("excitation_level", data=0) + file.create_dataset("temperature", data=T_KELVIN).attrs["unit"] = "K" + file.create_dataset("atomic_number", data=Z) + file.create_dataset("mass_number", data=A) + file.create_dataset("atomic_weight_ratio", data=awr) + file.create_dataset("radiation_length", data=radiation_length).attrs["unit"] = "g/cm2" + file.create_dataset("fissionable", data=False) + + written = write_stopping_power(file, pstar_dir, symbol, verbose=verbose) + file.close() + + if not written: + # No PSTAR data either — remove the empty file and signal failure + os.remove(out_path) + return None + + return mcdc_name + + +# -- Main ---------------------------------------------------------------------- + +def main(): + parser = argparse.ArgumentParser( + description="MC/DC proton data generator" + ) + parser.add_argument("--rewrite", action="store_true", default=False) + parser.add_argument("--verbose", action="store_true", default=False) + args = parser.parse_args() + rewrite = args.rewrite + verbose = args.verbose + + output_dir = os.getenv("MCDC_LIB_PROTON") + ace_dir = os.getenv("MCDC_ACELIB_PROTON") + pstar_dir = os.getenv("MCDC_PSTAR_LIB") + if ace_dir is None: + print_error("Environment variable $MCDC_ACELIB_PROTON is not set.") + if pstar_dir is None: + print_error("Environment variable $MCDC_PSTAR_LIB is not set.") + if output_dir is None: + print_error("Environment variable $MCDC_LIB_PROTON is not set.") + + os.makedirs(output_dir, exist_ok=True) + print(f"\nACE directory : {ace_dir}") + print(f"PSTAR directory : {pstar_dir}") + print(f"Output directory: {output_dir}\n") + + ace_files = sorted(f for f in os.listdir(ace_dir) if f.endswith(".ace")) + + # ── Pass 1: ACE files ───────────────────────────────────────────────────── + + if rewrite: + target_files = ace_files + else: + target_files = [] + for fname in ace_files: + try: + with open(os.path.join(ace_dir, fname)) as f: + hdr = ACEtk.Header.from_string(f.readline()) + Z, A, S, _ = decode_ace_zaid(hdr.zaid) + symbol = Z_TO_SYMBOL.get(Z, f"Z{Z}") + nuclide_name = f"{symbol}{A}" if S == 0 else f"{symbol}{A}m{S}" + if not any( + f.startswith(nuclide_name + "-") + for f in os.listdir(output_dir) + ): + target_files.append(fname) + except Exception: + target_files.append(fname) + + errors = [] + pbar = tqdm(target_files, disable=verbose, + bar_format="{l_bar}{bar}| {n_fmt}/{total_fmt} {postfix}") + + for ace_name in pbar: + pbar.set_postfix_str(ace_name) + try: + out = process_ace_file( + os.path.join(ace_dir, ace_name), + output_dir, + pstar_dir=pstar_dir, + verbose=verbose, + ) + if verbose: + print(f" -> wrote {out}") + except Exception as exc: + errors.append((ace_name, str(exc))) + if verbose: + import traceback + traceback.print_exc() + + # ── Pass 2: PSTAR-only isotopes (e.g. H, He) ───────────────────────────── + # For each entry in PSTAR_ONLY_ISOTOPES, create a stopping-power-only HDF5 + # file if one doesn't already exist (or if --rewrite is set). + + if pstar_dir is not None: + existing = set(os.listdir(output_dir)) + for symbol, A, awr in PSTAR_ONLY_ISOTOPES: + nuclide_name = f"{symbol}{A}" + mcdc_name = f"{nuclide_name}-{T_KELVIN}K.h5" + if not rewrite and mcdc_name in existing: + continue + try: + out = process_pstar_only_file( + symbol, A, awr, output_dir, pstar_dir, + verbose=verbose + ) + if out is None: + if verbose: + print(f" [skip] No PSTAR data for {nuclide_name}") + elif verbose: + print(f" -> wrote {out} [stopping power only]") + except Exception as exc: + errors.append((nuclide_name, str(exc))) + if verbose: + import traceback + traceback.print_exc() + + # ── Summary ─────────────────────────────────────────────────────────────── + + n_total = len(target_files) + len(PSTAR_ONLY_ISOTOPES) + print(f"\nDone. {n_total - len(errors)} succeeded, {len(errors)} failed.") + if errors: + print("\nFailed files:") + for name, msg in errors: + print(f" {name}: {msg}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/tools/data_library_generator/proton/water.py b/tools/data_library_generator/proton/water.py new file mode 100644 index 000000000..92799a0c4 --- /dev/null +++ b/tools/data_library_generator/proton/water.py @@ -0,0 +1,26 @@ +import h5py +import numpy as np + +def load_pstar_file(filepath): + energies, sps = [], [] + with open(filepath) as f: + for line in f: + parts = line.strip().split() + if len(parts) != 2: + continue + try: + energies.append(float(parts[0])) + sps.append(float(parts[1])) + except ValueError: + continue + return np.array(energies), np.array(sps) + +file = h5py.File("../../../proton_generated_lib/p_in_H2O.h5", "w") +pstar_path = "../../../pstar_lib/H2O.txt" +E_s, S_s = load_pstar_file(pstar_path) +X0 = 36.08 +sp = file.create_group("stopping_power") +sp.create_dataset("energy", data=E_s).attrs["unit"] = "MeV" +sp.create_dataset("total_stopping_power", data=S_s).attrs["unit"] = "MeV cm2/g" +rad_length = file.create_group("radiation_length") +rad_length.create_dataset("radiation_length", data=X0).attrs["unit"] = "g/cm2" \ No newline at end of file