diff --git a/mcdc/constant.py b/mcdc/constant.py index 5e3d02fa..a085c2d6 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -27,6 +27,7 @@ SCORE_SPACE_MOMENT_CURRENT = 10 SCORE_TIME_MOMENT_MU_SQ = 11 SCORE_SPACE_MOMENT_MU_SQ = 12 +SCORE_EDEP = 13 # Tally multipliers MULTIPLIER_ENERGY = 0 @@ -99,9 +100,17 @@ REACTION_NEUTRON_FISSION = 4 REACTION_NEUTRON_FISSION_PROMPT = 5 REACTION_NEUTRON_FISSION_DELAYED = 6 +# +REACTION_ELECTRON_ELASTIC_SCATTERING = 101 +REACTION_ELECTRON_ELASTIC_SMALL_ANGLE = 102 +REACTION_ELECTRON_ELASTIC_LARGE_ANGLE = 103 +REACTION_ELECTRON_IONIZATION = 104 +REACTION_ELECTRON_BREMSSTRAHLUNG = 105 +REACTION_ELECTRON_EXCITATION = 106 # Particle types PARTICLE_NEUTRON = 0 +PARTICLE_ELECTRON = 1 # Data DATA_NONE = 0 @@ -167,15 +176,34 @@ AXIS_Z = 2 AXIS_T = 3 +# ====================================================================================== # Physics +# ====================================================================================== + +# Constants LIGHT_SPEED = 2.99792458e10 # cm/s -NEUTRON_MASS = 939.565413e6 # eV/c^2 BOLTZMANN_K = 8.61733326e-5 # eV/K + +# Adjustables THERMAL_THRESHOLD_FACTOR = 400 +# Neutron properties +NEUTRON_MASS = 939.565413e6 # eV/c^2 + +# Electron energy +ELECTRON_MASS = 510.99895069e3 # eV/c^2 +ELECTRON_CUTOFF_ENERGY = 100 # eV (TODO: Is this adjustable?) +FINE_STRUCTURE_CONSTANT = 7.2973525693e-3 # ~1/137 +MU_CUTOFF = 0.999999 + +# ====================================================================================== +# Techniques +# ====================================================================================== + # Weight Windows Methods WW_USER = 0 WW_PREVIOUS = 1 + # Weight Windows Modifications WW_MIN = 0 WW_WOLLABER = 1 diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py index 31cde0a6..8e1682bf 100644 --- a/mcdc/mcdc_get/__init__.py +++ b/mcdc/mcdc_get/__init__.py @@ -44,6 +44,16 @@ import mcdc.mcdc_get.tabulated_energy_angle_distribution as tabulated_energy_angle_distribution +import mcdc.mcdc_get.element as element + +import mcdc.mcdc_get.electron_bremsstrahlung_reaction as electron_bremsstrahlung_reaction + +import mcdc.mcdc_get.electron_elastic_scattering_reaction as electron_elastic_scattering_reaction + +import mcdc.mcdc_get.electron_excitation_reaction as electron_excitation_reaction + +import mcdc.mcdc_get.electron_ionization_reaction as electron_ionization_reaction + import mcdc.mcdc_get.gpu_meta as gpu_meta import mcdc.mcdc_get.native_material as native_material diff --git a/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py b/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_bremsstrahlung_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_elastic_scattering_reaction.py b/mcdc/mcdc_get/electron_elastic_scattering_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_elastic_scattering_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_excitation_reaction.py b/mcdc/mcdc_get/electron_excitation_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_get/electron_excitation_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_get/electron_ionization_reaction.py b/mcdc/mcdc_get/electron_ionization_reaction.py new file mode 100644 index 00000000..d05a69d9 --- /dev/null +++ b/mcdc/mcdc_get/electron_ionization_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def binding_energy(index, electron_ionization_reaction, data): + offset = electron_ionization_reaction["binding_energy_offset"] + return data[offset + index] + + +@njit +def binding_energy_all(electron_ionization_reaction, data): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + return data[start:end] + + +@njit +def binding_energy_last(electron_ionization_reaction, data): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + return data[end - 1] + + +@njit +def binding_energy_chunk(start, length, electron_ionization_reaction, data): + start += electron_ionization_reaction["binding_energy_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/element.py b/mcdc/mcdc_get/element.py new file mode 100644 index 00000000..33f86e4f --- /dev/null +++ b/mcdc/mcdc_get/element.py @@ -0,0 +1,322 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs_energy_grid(index, element, data): + offset = element["xs_energy_grid_offset"] + return data[offset + index] + + +@njit +def xs_energy_grid_all(element, data): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + return data[start:end] + + +@njit +def xs_energy_grid_last(element, data): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + return data[end - 1] + + +@njit +def xs_energy_grid_chunk(start, length, element, data): + start += element["xs_energy_grid_offset"] + end = start + length + return data[start:end] + + +@njit +def total_xs(index, element, data): + offset = element["total_xs_offset"] + return data[offset + index] + + +@njit +def total_xs_all(element, data): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + return data[start:end] + + +@njit +def total_xs_last(element, data): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def total_xs_chunk(start, length, element, data): + start += element["total_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_xs(index, element, data): + offset = element["ionization_xs_offset"] + return data[offset + index] + + +@njit +def ionization_xs_all(element, data): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + return data[start:end] + + +@njit +def ionization_xs_last(element, data): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def ionization_xs_chunk(start, length, element, data): + start += element["ionization_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def elastic_xs(index, element, data): + offset = element["elastic_xs_offset"] + return data[offset + index] + + +@njit +def elastic_xs_all(element, data): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + return data[start:end] + + +@njit +def elastic_xs_last(element, data): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def elastic_xs_chunk(start, length, element, data): + start += element["elastic_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def excitation_xs(index, element, data): + offset = element["excitation_xs_offset"] + return data[offset + index] + + +@njit +def excitation_xs_all(element, data): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + return data[start:end] + + +@njit +def excitation_xs_last(element, data): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def excitation_xs_chunk(start, length, element, data): + start += element["excitation_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def bremsstrahlung_xs(index, element, data): + offset = element["bremsstrahlung_xs_offset"] + return data[offset + index] + + +@njit +def bremsstrahlung_xs_all(element, data): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + return data[start:end] + + +@njit +def bremsstrahlung_xs_last(element, data): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + return data[end - 1] + + +@njit +def bremsstrahlung_xs_chunk(start, length, element, data): + start += element["bremsstrahlung_xs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_reaction_IDs(index, element, data): + offset = element["ionization_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def ionization_reaction_IDs_all(element, data): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + return data[start:end] + + +@njit +def ionization_reaction_IDs_last(element, data): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + return data[end - 1] + + +@njit +def ionization_reaction_IDs_chunk(start, length, element, data): + start += element["ionization_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def elastic_scattering_reaction_IDs(index, element, data): + offset = element["elastic_scattering_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def elastic_scattering_reaction_IDs_all(element, data): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + return data[start:end] + + +@njit +def elastic_scattering_reaction_IDs_last(element, data): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + return data[end - 1] + + +@njit +def elastic_scattering_reaction_IDs_chunk(start, length, element, data): + start += element["elastic_scattering_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def excitation_reaction_IDs(index, element, data): + offset = element["excitation_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def excitation_reaction_IDs_all(element, data): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + return data[start:end] + + +@njit +def excitation_reaction_IDs_last(element, data): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + return data[end - 1] + + +@njit +def excitation_reaction_IDs_chunk(start, length, element, data): + start += element["excitation_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def bremsstrahlung_reaction_IDs(index, element, data): + offset = element["bremsstrahlung_reaction_IDs_offset"] + return data[offset + index] + + +@njit +def bremsstrahlung_reaction_IDs_all(element, data): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + return data[start:end] + + +@njit +def bremsstrahlung_reaction_IDs_last(element, data): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + return data[end - 1] + + +@njit +def bremsstrahlung_reaction_IDs_chunk(start, length, element, data): + start += element["bremsstrahlung_reaction_IDs_offset"] + end = start + length + return data[start:end] + + +@njit +def ionization_subshell_binding_energy(index, element, data): + offset = element["ionization_subshell_binding_energy_offset"] + return data[offset + index] + + +@njit +def ionization_subshell_binding_energy_all(element, data): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + return data[start:end] + + +@njit +def ionization_subshell_binding_energy_last(element, data): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + return data[end - 1] + + +@njit +def ionization_subshell_binding_energy_chunk(start, length, element, data): + start += element["ionization_subshell_binding_energy_offset"] + end = start + length + return data[start:end] diff --git a/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py b/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py index 1677db2b..5cf2f0dd 100644 --- a/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py +++ b/mcdc/mcdc_get/neutron_inelastic_scattering_reaction.py @@ -26,9 +26,7 @@ def spectrum_probability_grid_last(neutron_inelastic_scattering_reaction, data): @njit -def spectrum_probability_grid_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def spectrum_probability_grid_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] end = start + length return data[start:end] @@ -51,9 +49,7 @@ def spectrum_probability(index_1, index_2, neutron_inelastic_scattering_reaction @njit -def spectrum_probability_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def spectrum_probability_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["spectrum_probability_offset"] end = start + length return data[start:end] @@ -82,9 +78,7 @@ def energy_spectrum_IDs_last(neutron_inelastic_scattering_reaction, data): @njit -def energy_spectrum_IDs_chunk( - start, length, neutron_inelastic_scattering_reaction, data -): +def energy_spectrum_IDs_chunk(start, length, neutron_inelastic_scattering_reaction, data): start += neutron_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] end = start + length return data[start:end] diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py index 4df13dfc..5664c21f 100644 --- a/mcdc/mcdc_set/__init__.py +++ b/mcdc/mcdc_set/__init__.py @@ -44,6 +44,16 @@ import mcdc.mcdc_set.tabulated_energy_angle_distribution as tabulated_energy_angle_distribution +import mcdc.mcdc_set.element as element + +import mcdc.mcdc_set.electron_bremsstrahlung_reaction as electron_bremsstrahlung_reaction + +import mcdc.mcdc_set.electron_elastic_scattering_reaction as electron_elastic_scattering_reaction + +import mcdc.mcdc_set.electron_excitation_reaction as electron_excitation_reaction + +import mcdc.mcdc_set.electron_ionization_reaction as electron_ionization_reaction + import mcdc.mcdc_set.gpu_meta as gpu_meta import mcdc.mcdc_set.native_material as native_material diff --git a/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py b/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_bremsstrahlung_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_elastic_scattering_reaction.py b/mcdc/mcdc_set/electron_elastic_scattering_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_elastic_scattering_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_excitation_reaction.py b/mcdc/mcdc_set/electron_excitation_reaction.py new file mode 100644 index 00000000..fdbf8e75 --- /dev/null +++ b/mcdc/mcdc_set/electron_excitation_reaction.py @@ -0,0 +1,3 @@ +# The following is automatically generated by code_factory.py + +from numba import njit diff --git a/mcdc/mcdc_set/electron_ionization_reaction.py b/mcdc/mcdc_set/electron_ionization_reaction.py new file mode 100644 index 00000000..df76b1b9 --- /dev/null +++ b/mcdc/mcdc_set/electron_ionization_reaction.py @@ -0,0 +1,32 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def binding_energy(index, electron_ionization_reaction, data, value): + offset = electron_ionization_reaction["binding_energy_offset"] + data[offset + index] = value + + +@njit +def binding_energy_all(electron_ionization_reaction, data, value): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + data[start:end] = value + + +@njit +def binding_energy_last(electron_ionization_reaction, data, value): + start = electron_ionization_reaction["binding_energy_offset"] + size = electron_ionization_reaction["binding_energy_length"] + end = start + size + data[end - 1] = value + + +@njit +def binding_energy_chunk(start, length, electron_ionization_reaction, data, value): + start += electron_ionization_reaction["binding_energy_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/element.py b/mcdc/mcdc_set/element.py new file mode 100644 index 00000000..11a248a9 --- /dev/null +++ b/mcdc/mcdc_set/element.py @@ -0,0 +1,322 @@ +# The following is automatically generated by code_factory.py + +from numba import njit + + +@njit +def xs_energy_grid(index, element, data, value): + offset = element["xs_energy_grid_offset"] + data[offset + index] = value + + +@njit +def xs_energy_grid_all(element, data, value): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + data[start:end] = value + + +@njit +def xs_energy_grid_last(element, data, value): + start = element["xs_energy_grid_offset"] + size = element["xs_energy_grid_length"] + end = start + size + data[end - 1] = value + + +@njit +def xs_energy_grid_chunk(start, length, element, data, value): + start += element["xs_energy_grid_offset"] + end = start + length + data[start:end] = value + + +@njit +def total_xs(index, element, data, value): + offset = element["total_xs_offset"] + data[offset + index] = value + + +@njit +def total_xs_all(element, data, value): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def total_xs_last(element, data, value): + start = element["total_xs_offset"] + size = element["total_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def total_xs_chunk(start, length, element, data, value): + start += element["total_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_xs(index, element, data, value): + offset = element["ionization_xs_offset"] + data[offset + index] = value + + +@njit +def ionization_xs_all(element, data, value): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def ionization_xs_last(element, data, value): + start = element["ionization_xs_offset"] + size = element["ionization_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_xs_chunk(start, length, element, data, value): + start += element["ionization_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def elastic_xs(index, element, data, value): + offset = element["elastic_xs_offset"] + data[offset + index] = value + + +@njit +def elastic_xs_all(element, data, value): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def elastic_xs_last(element, data, value): + start = element["elastic_xs_offset"] + size = element["elastic_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def elastic_xs_chunk(start, length, element, data, value): + start += element["elastic_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def excitation_xs(index, element, data, value): + offset = element["excitation_xs_offset"] + data[offset + index] = value + + +@njit +def excitation_xs_all(element, data, value): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def excitation_xs_last(element, data, value): + start = element["excitation_xs_offset"] + size = element["excitation_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def excitation_xs_chunk(start, length, element, data, value): + start += element["excitation_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def bremsstrahlung_xs(index, element, data, value): + offset = element["bremsstrahlung_xs_offset"] + data[offset + index] = value + + +@njit +def bremsstrahlung_xs_all(element, data, value): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + data[start:end] = value + + +@njit +def bremsstrahlung_xs_last(element, data, value): + start = element["bremsstrahlung_xs_offset"] + size = element["bremsstrahlung_xs_length"] + end = start + size + data[end - 1] = value + + +@njit +def bremsstrahlung_xs_chunk(start, length, element, data, value): + start += element["bremsstrahlung_xs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_reaction_IDs(index, element, data, value): + offset = element["ionization_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def ionization_reaction_IDs_all(element, data, value): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + data[start:end] = value + + +@njit +def ionization_reaction_IDs_last(element, data, value): + start = element["ionization_reaction_IDs_offset"] + size = element["N_ionization_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_reaction_IDs_chunk(start, length, element, data, value): + start += element["ionization_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def elastic_scattering_reaction_IDs(index, element, data, value): + offset = element["elastic_scattering_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def elastic_scattering_reaction_IDs_all(element, data, value): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + data[start:end] = value + + +@njit +def elastic_scattering_reaction_IDs_last(element, data, value): + start = element["elastic_scattering_reaction_IDs_offset"] + size = element["N_elastic_scattering_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def elastic_scattering_reaction_IDs_chunk(start, length, element, data, value): + start += element["elastic_scattering_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def excitation_reaction_IDs(index, element, data, value): + offset = element["excitation_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def excitation_reaction_IDs_all(element, data, value): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + data[start:end] = value + + +@njit +def excitation_reaction_IDs_last(element, data, value): + start = element["excitation_reaction_IDs_offset"] + size = element["N_excitation_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def excitation_reaction_IDs_chunk(start, length, element, data, value): + start += element["excitation_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def bremsstrahlung_reaction_IDs(index, element, data, value): + offset = element["bremsstrahlung_reaction_IDs_offset"] + data[offset + index] = value + + +@njit +def bremsstrahlung_reaction_IDs_all(element, data, value): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + data[start:end] = value + + +@njit +def bremsstrahlung_reaction_IDs_last(element, data, value): + start = element["bremsstrahlung_reaction_IDs_offset"] + size = element["N_bremsstrahlung_reaction"] + end = start + size + data[end - 1] = value + + +@njit +def bremsstrahlung_reaction_IDs_chunk(start, length, element, data, value): + start += element["bremsstrahlung_reaction_IDs_offset"] + end = start + length + data[start:end] = value + + +@njit +def ionization_subshell_binding_energy(index, element, data, value): + offset = element["ionization_subshell_binding_energy_offset"] + data[offset + index] = value + + +@njit +def ionization_subshell_binding_energy_all(element, data, value): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + data[start:end] = value + + +@njit +def ionization_subshell_binding_energy_last(element, data, value): + start = element["ionization_subshell_binding_energy_offset"] + size = element["ionization_subshell_binding_energy_length"] + end = start + size + data[end - 1] = value + + +@njit +def ionization_subshell_binding_energy_chunk(start, length, element, data, value): + start += element["ionization_subshell_binding_energy_offset"] + end = start + length + data[start:end] = value diff --git a/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py b/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py index 65509585..81468be2 100644 --- a/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py +++ b/mcdc/mcdc_set/neutron_inelastic_scattering_reaction.py @@ -4,9 +4,7 @@ @njit -def spectrum_probability_grid( - index, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_grid(index, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] data[offset + index] = value @@ -28,18 +26,14 @@ def spectrum_probability_grid_last(neutron_inelastic_scattering_reaction, data, @njit -def spectrum_probability_grid_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_grid_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["spectrum_probability_grid_offset"] end = start + length data[start:end] = value @njit -def spectrum_probability_vector( - index_1, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_vector(index_1, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_offset"] stride = neutron_inelastic_scattering_reaction["N_spectrum"] start = offset + index_1 * stride @@ -48,18 +42,14 @@ def spectrum_probability_vector( @njit -def spectrum_probability( - index_1, index_2, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability(index_1, index_2, neutron_inelastic_scattering_reaction, data, value): offset = neutron_inelastic_scattering_reaction["spectrum_probability_offset"] stride = neutron_inelastic_scattering_reaction["N_spectrum"] data[offset + index_1 * stride + index_2] = value @njit -def spectrum_probability_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def spectrum_probability_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["spectrum_probability_offset"] end = start + length data[start:end] = value @@ -88,9 +78,7 @@ def energy_spectrum_IDs_last(neutron_inelastic_scattering_reaction, data, value) @njit -def energy_spectrum_IDs_chunk( - start, length, neutron_inelastic_scattering_reaction, data, value -): +def energy_spectrum_IDs_chunk(start, length, neutron_inelastic_scattering_reaction, data, value): start += neutron_inelastic_scattering_reaction["energy_spectrum_IDs_offset"] end = start + length data[start:end] = value diff --git a/mcdc/mcdc_set/tabulated_energy_angle_distribution.py b/mcdc/mcdc_set/tabulated_energy_angle_distribution.py index 97040280..e25747d1 100644 --- a/mcdc/mcdc_set/tabulated_energy_angle_distribution.py +++ b/mcdc/mcdc_set/tabulated_energy_angle_distribution.py @@ -171,9 +171,7 @@ def cosine_offset__last(tabulated_energy_angle_distribution, data, value): @njit -def cosine_offset__chunk( - start, length, tabulated_energy_angle_distribution, data, value -): +def cosine_offset__chunk(start, length, tabulated_energy_angle_distribution, data, value): start += tabulated_energy_angle_distribution["cosine_offset__offset"] end = start + length data[start:end] = value diff --git a/mcdc/object_/element.py b/mcdc/object_/element.py new file mode 100644 index 00000000..5e808aab --- /dev/null +++ b/mcdc/object_/element.py @@ -0,0 +1,178 @@ +import h5py +import numpy as np +import os + +from numpy import float64 +from numpy.typing import NDArray + +#### + +from mcdc.object_.base import ObjectNonSingleton +from mcdc.object_.reaction import ( + ReactionElectronElasticScattering, + ReactionElectronExcitation, + ReactionElectronBremsstrahlung, + ReactionElectronIonization, +) +from mcdc.print_ import print_1d_array, print_error + + +# ====================================================================================== +# Element +# ====================================================================================== + + +class Element(ObjectNonSingleton): + # Annotations for Numba mode + label: str = "element" + # + name: str + atomic_weight_ratio: float + atomic_number: int + # XS on xs_energy_grid + xs_energy_grid: NDArray[float64] + total_xs: NDArray[float64] + ionization_xs: NDArray[float64] + elastic_xs: NDArray[float64] + excitation_xs: NDArray[float64] + bremsstrahlung_xs: NDArray[float64] + # Reactions + ionization_reactions: list[ReactionElectronIonization] + elastic_scattering_reactions: list[ReactionElectronElasticScattering] + excitation_reactions: list[ReactionElectronExcitation] + bremsstrahlung_reactions: list[ReactionElectronBremsstrahlung] + # Ionization element attributes + ionization_subshell_binding_energy: NDArray[float64] + + def __init__(self, element_name: str): + super().__init__() + + self.name = element_name + + # Set attributes from the hdf5 file + dir_name = os.getenv("MCDC_LIB") + file_name = f"{element_name}.h5" + file = h5py.File(f"{dir_name}/{file_name}", "r") + + # Basic properties + self.atomic_weight_ratio = float(file["atomic_weight_ratio"][()]) + self.atomic_number = int(file["atomic_number"][()]) + + # The reactions + rx_names = [ + "elastic_scattering", + "excitation", + "bremsstrahlung", + "ionization", + ] + + # The reaction MTs + MTs = {} + for name in rx_names: + if name not in file["electron_reactions"]: + MTs[name] = [] + continue + + MTs[name] = [ + x for x in file[f"electron_reactions/{name}"] if x.startswith("MT") + ] + + # ========================================================================== + # Reaction XS + # ========================================================================== + + # Energy grid + xs_energy = file["electron_reactions/xs_energy_grid"][()] + self.xs_energy_grid = xs_energy + + # The total XS + self.total_xs = np.zeros_like(self.xs_energy_grid) + self.elastic_xs = np.zeros_like(self.xs_energy_grid) + self.excitation_xs = np.zeros_like(self.xs_energy_grid) + self.bremsstrahlung_xs = np.zeros_like(self.xs_energy_grid) + self.ionization_xs = np.zeros_like(self.xs_energy_grid) + + xs_containers = [ + self.elastic_xs, + self.excitation_xs, + self.bremsstrahlung_xs, + self.ionization_xs, + ] + for xs_container, rx_name in list(zip(xs_containers, rx_names)): + for MT in MTs[rx_name]: + xs = file[f"electron_reactions/{rx_name}/{MT}/xs"] + xs_container[xs.attrs["offset"] :] += xs[()] + + # Total = sum of xs' + self.total_xs = ( + self.elastic_xs + + self.excitation_xs + + self.bremsstrahlung_xs + + self.ionization_xs + ) + + # ========================================================================== + # The reactions + # ========================================================================== + + self.elastic_scattering_reactions = [] + self.excitation_reactions = [] + self.bremsstrahlung_reactions = [] + self.ionization_reactions = [] + + rx_containers = [ + self.elastic_scattering_reactions, + self.excitation_reactions, + self.bremsstrahlung_reactions, + self.ionization_reactions, + ] + rx_classes = [ + ReactionElectronElasticScattering, + ReactionElectronExcitation, + ReactionElectronBremsstrahlung, + ReactionElectronIonization, + ] + + for rx_container, rx_name, rx_class in list(zip(rx_containers, rx_names, rx_classes)): + for MT in MTs[rx_name]: + h5_group = file[f"electron_reactions/{rx_name}/{MT}"] + reaction = rx_class.from_h5_group(h5_group) + rx_container.append(reaction) + + # ========================================================================== + # Ionization element attributes + # ========================================================================== + + be_list = [] + if len(MTs["ionization"]) > 0: + g = file[f"electron_reactions/ionization/{MTs['ionization'][0]}"] + subshells = g["subshells"] + subshell_names = list(subshells.keys()) + for name in subshell_names: + sh = subshells[name] + be_list.append(float(sh["binding_energy"][()])) + + self.ionization_subshell_binding_energy = np.asarray(be_list) + + file.close() + + def __repr__(self): + text = "\n" + text += f"Element\n" + text += f" - ID: {self.ID}\n" + text += f" - Name: {self.name}\n" + text += f" - Atomic number: {self.atomic_number}\n" + text += f" - Atomic weight ratio: {self.atomic_weight_ratio}\n" + text += f" - Reaction MTs\n" + text += f" - Ionization: {[int(x.MT) for x in self.ionization_reactions]}\n" + text += f" - Elastic scattering: {[int(x.MT) for x in self.elastic_scattering_reactions]}\n" + text += f" - Excitation: {[int(x.MT) for x in self.excitation_reactions]}\n" + text += f" - Bremsstrahlung: {[int(x.MT) for x in self.bremsstrahlung_reactions]}\n" + text += f" - Reaction cross-sections (eV, barns)\n" + text += f" - Energy grid {print_1d_array(self.xs_energy_grid)}\n" + text += f" - Total {print_1d_array(self.total_xs)}\n" + text += f" - Ionization {print_1d_array(self.ionization_xs)}\n" + text += f" - Elastic scattering {print_1d_array(self.elastic_xs)}\n" + text += f" - Excitation {print_1d_array(self.excitation_xs)}\n" + text += f" - Bremsstrahlung {print_1d_array(self.bremsstrahlung_xs)}\n" + return text \ No newline at end of file diff --git a/mcdc/object_/material.py b/mcdc/object_/material.py index bf3c4831..b8549007 100644 --- a/mcdc/object_/material.py +++ b/mcdc/object_/material.py @@ -11,8 +11,10 @@ from mcdc.constant import MATERIAL, MATERIAL_MG from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.nuclide import Nuclide +from mcdc.object_.element import Element from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error +from mcdc.object_.util import ISOTOPIC_ABUNDANCE # ====================================================================================== # Material base class @@ -89,17 +91,22 @@ class Material(MaterialBase): # Annotations for Numba mode label: str = "native_material" non_numba: list[str] = ["nuclide_composition"] + non_numba_element: list[str] = ["element_composition"] # nuclide_composition: dict[Nuclide, float] + element_composition: dict[Element, float] # nuclides: list[Nuclide] + elements: list[Element] nuclide_densities: NDArray[float64] + element_densities: NDArray[float64] def __init__( self, name: str = "", nuclide_composition: dict[str, float] = {}, - temperature: float = 293.6, + element_composition: dict[str, float] = {}, + temperature: float = 293.6 ): type_ = MATERIAL super().__init__(type_, name) @@ -110,15 +117,73 @@ def __init__( # Dictionary connecting nuclides to respective densities self.nuclide_composition = {} + # Dictionary connecting elements to respective densities + self.element_composition = {} + # Numba representation of nuclide_composition self.nuclides = [] self.nuclide_densities = np.zeros(len(nuclide_composition)) + # Numba representation of element_composition + self.elements = [] + self.element_densities = np.zeros(len(element_composition)) + + # ============================================================================== + # Argument checks + # ============================================================================== + + # Check that either nuclide or element composition is provided, not both + if len(nuclide_composition) > 0 and len(element_composition) > 0: + print_error("Cannot specify both nuclide_composition and element_composition") + + if len(nuclide_composition) == 0 and len(element_composition) == 0: + print_error("Must specify either nuclide_composition or element_composition") + # Check if library directory is set lib_dir = os.getenv("MCDC_LIB") if lib_dir is None: print_error("Environment variable MCDC_LIB is not set") + # ============================================================================== + # Element-based construction + # ============================================================================== + + # Loop over the items in the composition + for i, (key, value) in enumerate(element_composition.items()): + element_name = key + element_density = value + temperature = 0.0 # Elements are always at 0 K + + # Check if element is available in the library + file_name = f"{element_name}.h5" + if not file_name in os.listdir(lib_dir): + print_error( + f"Element {element_name} is not available in the library" + ) + + # Check if element is already created + found = False + for element in simulation.elements: + if ( + element.name == element_name + ): + found = True + break + + # Create the element to objects if needed + if not found: + element = Element(element_name) + + # Register the element composition + self.elements.append(element) + self.element_densities[i] = element_density + self.element_composition[element] = element_density + + + # ============================================================================== + # Nuclide-based construction + # ============================================================================== + # Loop over the items in the composition for i, (key, value) in enumerate(nuclide_composition.items()): nuclide_name = key @@ -157,17 +222,27 @@ def __init__( if nuclide.fissionable: self.fissionable = True + # Create undefined compositions + if self.nuclides == []: + set_nuclides_from_elements(self) + if self.elements == []: + set_elements_from_nuclides(self) + def __repr__(self): text = super().__repr__() text += f" - Temperature: {self.temperature} K\n" text += f" - Nuclide composition [atoms/barn-cm]\n" for nuclide in self.nuclide_composition.keys(): text += f" - {nuclide.name:<5} | {self.nuclide_composition[nuclide]}\n" + text += f" - Element composition [atoms/barn-cm]\n" + for element in self.element_composition.keys(): + text += f" - {element.name:<5} | {self.element_composition[element]}\n" return text # Currently supported temperatures TEMPERATURES = [0.1, 233.15, 273.15, 293.6, 600.0, 900.0, 1200.0, 2500.0] +TEMPERATURES_ELEMENT = [0.0] # ====================================================================================== @@ -389,3 +464,110 @@ def __repr__(self): text += f" - speed {print_1d_array(self.mgxs_speed)}\n" text += f" - lambda {print_1d_array(self.mgxs_decay_rate)}\n" return text + + +# ====================================================================================== +# Helpers +# ====================================================================================== + +def set_nuclides_from_elements(material): + material.nuclides = [] + material.nuclide_densities = [] + material.nuclide_composition = {} + + # Check if library directory is set + lib_dir = os.getenv("MCDC_LIB") + if lib_dir is None: + print_error("Environment variable MCDC_LIB is not set") + + # Get supported temperature + nearest_temperature = min( + TEMPERATURES, key=lambda x: abs(x - material.temperature) + ) + + for element, element_density in material.element_composition.items(): + # To make sure that the abundance is normalized + norm = 0.0 + for (_, abundance) in ISOTOPIC_ABUNDANCE[element.name].items(): + norm += abundance + + # Loop over the nuclide composition + for nuclide_name, abundance in ISOTOPIC_ABUNDANCE[element.name].items(): + # Check if nuclide-temperature is available in the library + file_name = f"{nuclide_name}-{nearest_temperature}K.h5" + if not file_name in os.listdir(lib_dir): + print_error( + f"Nuclide {nuclide_name} at temperature {nearest_temperature} K is not available in the library" + ) + + # Check if nuclide is already created + found = False + for nuclide in simulation.nuclides: + if ( + nuclide.name == nuclide_name + and nearest_temperature == nuclide.temperature + ): + found = True + break + + # Create the nuclide to objects if needed + if not found: + nuclide = Nuclide(nuclide_name, nearest_temperature) + + # Calculate nuclide density + nuclide_density = element_density * abundance / norm + + # Register the nuclide composition + material.nuclides.append(nuclide) + material.nuclide_densities.append(nuclide_density) + material.nuclide_composition[nuclide] = nuclide_density + + # Some flags + if nuclide.fissionable: + material.fissionable = True + + material.nuclide_densities = np.array(material.nuclide_densities) + + +def set_elements_from_nuclides(material): + material.elements = [] + material.element_densities = [] + material.element_composition = {} + + # Get the list of the element names + element_names = [] + for nuclide in material.nuclides: + name = nuclide.name[:2] + if name[1].isdigit(): + name = name[0] + if name not in element_names: + element_names.append(name) + material.element_densities = np.zeros(len(element_names)) + + # Iterate over all named elements + for i, element_name in enumerate(element_names): + # Check if element is already created + found = False + for element in simulation.elements: + if element.name == element_name: + found = True + break + + # Create the element to objects if needed + if not found: + element = Element(element_name) + + material.elements.append(element) + + # Iterate over all nuclides to get the total density + density = 0.0 + for nuclide, nuclide_density in material.nuclide_composition.items(): + # Skip if non-isotope + if nuclide.name[: len(element_name)] != element_name: + continue + + # Accumulate density + density += material.nuclide_composition[nuclide] + + material.element_densities[i] = density + material.element_composition[element] = density diff --git a/mcdc/object_/reaction.py b/mcdc/object_/reaction.py index 9e07710b..24c87675 100644 --- a/mcdc/object_/reaction.py +++ b/mcdc/object_/reaction.py @@ -4,6 +4,7 @@ #### +from mcdc.object_.data import DataBase, DataTable import mcdc.object_.distribution as distribution from mcdc.constant import ( @@ -16,8 +17,13 @@ REACTION_NEUTRON_ELASTIC_SCATTERING, REACTION_NEUTRON_FISSION, REACTION_NEUTRON_INELASTIC_SCATTERING, + REACTION_ELECTRON_BREMSSTRAHLUNG, + REACTION_ELECTRON_EXCITATION, + REACTION_ELECTRON_ELASTIC_SCATTERING, + REACTION_ELECTRON_IONIZATION, REFERENCE_FRAME_COM, REFERENCE_FRAME_LAB, + MU_CUTOFF, ) from mcdc.object_.base import ObjectPolymorphic from mcdc.object_.distribution import ( @@ -44,7 +50,7 @@ class ReactionBase(ObjectPolymorphic): # MT: int xs: NDArray[float64] - xs_offset_: int # "xs_offset" ir reserved for "xs" + xs_offset_: int # "xs_offset" is reserved for "xs" reference_frame: int def __init__(self, type_, MT, xs, xs_offset, reference_frame): @@ -73,6 +79,14 @@ def decode_type(type_): return "Neutron inelastic scattering" elif type_ == REACTION_NEUTRON_FISSION: return "Neutron fission" + elif type_ == REACTION_ELECTRON_IONIZATION: + return "Electron ionization" + elif type_ == REACTION_ELECTRON_ELASTIC_SCATTERING: + return "Electron elastic scattering" + elif type_ == REACTION_ELECTRON_BREMSSTRAHLUNG: + return "Electron bremsstrahlung" + elif type_ == REACTION_ELECTRON_EXCITATION: + return "Electron excitation" def decode_reference_frame(type_): @@ -282,6 +296,219 @@ def __repr__(self): return text +# ====================================================================================== +# Electron ionization +# ====================================================================================== + + +class ReactionElectronIonization(ReactionBase): + label: str = "electron_ionization_reaction" + # + N_subshell: int + subshell_xs: DistributionBase + binding_energy: NDArray[float64] + subshell_product: DistributionBase + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + subshell_xs, + binding_energy, + subshell_product, + ): + type_ = REACTION_ELECTRON_IONIZATION + super().__init__(type_, MT, xs, xs_offset, reference_frame) + + self.N_subshell = len(subshell_xs) + self.subshell_xs = subshell_xs + self.binding_energy = binding_energy + self.subshell_product = subshell_product + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + subshells = h5_group["subshells"] + subshell_names = list(subshells.keys()) + + subshell_xs = [] + subshell_product = [] + binding_energy = [] + + for name in subshell_names: + subshell = subshells[name] + + # Subshell cross section table (each has its own energy grid) + xs_sub = DataTable( + subshell["energy_grid"][()], + subshell["xs"][()] + ) + subshell_xs.append(xs_sub) + binding_energy.append(subshell[0].x[0]) + + # Secondary electron energy distribution + product = subshell["product"] + prod_dist = DistributionMultiTable( + product["energy_grid"][()], + product["energy_offset"][()], + product["value"][()], + product["PDF"][()] + ) + subshell_product.append(prod_dist) + + return cls( + MT, + xs, + xs_offset, + reference_frame, + subshell_xs, + binding_energy, + subshell_product, + ) + + def __repr__(self): + text = super().__repr__() + text += f" - Number of subshells: {self.N_subshell}\n" + for i in range(self.N_subshell): + text += f" - Subshell {i+1}\n" + text += f" - XS: DataTable [ID: {self.subshell_xs[i].ID}]\n" + prod = self.subshell_product[i] + text += ( + f" - Secondary electron spectrum: " + f"{distribution.decode_type(prod.type)} [ID: {prod.ID}]\n" + ) + return text + + +# ====================================================================================== +# Electron elastic scattering +# ====================================================================================== + + +class ReactionElectronElasticScattering(ReactionBase): + label: str = "electron_elastic_scattering_reaction" + # + xs_large: DataBase + mu: DistributionMultiTable + + def __init__( + self, + MT, + xs, + xs_offset, + reference_frame, + xs_large, + mu, + ): + type_ = REACTION_ELECTRON_ELASTIC_SCATTERING + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.mu_cut = MU_CUTOFF + self.xs_large = xs_large + self.mu = mu + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + # Large angle data + large_angle = h5_group["large_angle"] + + xs_large = DataTable( + large_angle["xs_energy"][()], + large_angle["xs"][()] + ) + + mu_g = large_angle["scattering_cosine"] + mu = DistributionMultiTable( + mu_g["energy_grid"][()], + mu_g["energy_offset"][()], + mu_g["value"][()], + mu_g["PDF"][()] + ) + + return cls( + MT, + xs, + xs_offset, + reference_frame, + xs_large, + mu, + ) + + def __repr__(self): + text = super().__repr__() + text += f" - Mu cut: {self.mu_cut}\n" + text += f" - Large angle XS: DataTable [ID: {self.xs_large.ID}]\n" + text += ( + f" - Scattering cosine: {distribution.decode_type(self.mu.type)} " + f"[ID: {self.mu.ID}]\n" + ) + return text + + +# ====================================================================================== +# Electron bremsstrahlung +# ====================================================================================== + + +class ReactionElectronBremsstrahlung(ReactionBase): + label: str = "electron_bremsstrahlung_reaction" + # + eloss: DataBase + + def __init__(self, MT, xs, xs_offset, reference_frame, eloss): + type_ = REACTION_ELECTRON_BREMSSTRAHLUNG + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.eloss = eloss + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + base = h5_group["energy_loss"] + eloss = DataTable(base["energy"][()], base["value"][()]) + + return cls(MT, xs, xs_offset, reference_frame, eloss) + + def __repr__(self): + text = super().__repr__() + text += f" - Energy loss: DataTable [ID: {self.eloss.ID}]\n" + return text + + +# ====================================================================================== +# Electron excitation +# ====================================================================================== + + +class ReactionElectronExcitation(ReactionBase): + label: str = "electron_excitation_reaction" + # + eloss: DataBase + + def __init__(self, MT, xs, xs_offset, reference_frame, eloss): + type_ = REACTION_ELECTRON_EXCITATION + super().__init__(type_, MT, xs, xs_offset, reference_frame) + self.eloss = eloss + + @classmethod + def from_h5_group(cls, h5_group): + MT, xs, xs_offset, reference_frame = set_basic_properties(h5_group) + + base = h5_group["energy_loss"] + eloss = DataTable(base["energy"][()], base["value"][()]) + + return cls(MT, xs, xs_offset, reference_frame, eloss) + + def __repr__(self): + text = super().__repr__() + text += f" - Energy loss: DataTable [ID: {self.eloss.ID}]\n" + return text + + # ====================================================================================== # Helper functions # ====================================================================================== @@ -395,4 +622,4 @@ def set_energy_distribution(h5_group): else: print_error(f"Unsupported energy spectrum of type {spectrum_type}") - return energy_spectrum + return energy_spectrum \ No newline at end of file diff --git a/mcdc/object_/simulation.py b/mcdc/object_/simulation.py index 95897f42..7048db5d 100644 --- a/mcdc/object_/simulation.py +++ b/mcdc/object_/simulation.py @@ -10,6 +10,7 @@ if TYPE_CHECKING: from mcdc.object_.cell import Cell, Region + from mcdc.object_.element import Element from mcdc.object_.material import MaterialBase from mcdc.object_.nuclide import Nuclide from mcdc.object_.reaction import ReactionBase @@ -56,6 +57,7 @@ class Simulation(ObjectSingleton): data: list[DataBase] distributions: list[DistributionBase] materials: list[MaterialBase] + elements: list[Element] nuclides: list[Nuclide] reactions: list[ReactionBase] sources: list[Source] @@ -138,6 +140,7 @@ def __init__(self): self.data = [DataNone()] self.distributions = [DistributionNone()] self.materials = [] + self.elements = [] self.nuclides = [] self.reactions = [] self.sources = [] diff --git a/mcdc/object_/util.py b/mcdc/object_/util.py index 2763a3f5..5f114c17 100644 --- a/mcdc/object_/util.py +++ b/mcdc/object_/util.py @@ -63,6 +63,471 @@ def is_sorted(a): return np.all(a[:-1] <= a[1:]) +# ====================================================================================== +# Data +# ====================================================================================== + +# Natural isotopic abundance data from +# https://www.nndc.bnl.gov/walletcards/search.html + +ISOTOPIC_ABUNDANCE = { + "H": { + "H1": (99.972 + 99.999) / 2, + "H2": (0.001 + 0.028) / 2, + }, + "He": { + "He3": 0.0002, + "He4": 99.9998, + }, + "Li": { + "Li6": (1.9 + 7.8) / 2, + "Li7": (92.2 + 98.1) / 2, + }, + "Be": { + "Be9": 100.0, + }, + "B": { + "B10": (18.9 + 20.4) / 2, + "B11": (79.6 + 81.1) / 2, + }, + "C": { + "C12": (98.84 + 99.04) / 2, + "C13": (0.96 + 1.16) / 2, + }, + "N": { + "N14": (99.578 + 99.663) / 2, + "N15": (0.337 + 0.422) / 2, + }, + "O": { + "O16": (99.738 + 99.776) / 2, + "O17": (0.0367 + 0.04) / 2, + "O18": (0.187 + 0.222) / 2, + }, + "F": { + "F19": 100.0, + }, + "Ne": { + "Ne20": 90.48, + "Ne21": 0.27, + "Ne22": 9.25, + }, + "Na": { + "Na23": 100.0, + }, + "Mg": { + "Mg24": (78.88 + 79.05) / 2, + "Mg25": (9.988 + 10.034) / 2, + "Mg26": (10.96 + 11.09) / 2, + }, + "Al": { + "Al27": 100.0, + }, + "Si": { + "Si28": (92.191 + 92.318) / 2, + "Si29": (4.645 + 4.699) / 2, + "Si30": (3.037 + 3.11) / 2, + }, + "P": { + "P31": 100.0, + }, + "S": { + "S32": (94.41 + 95.29) / 2, + "S33": (0.729 + 0.797) / 2, + "S34": (3.96 + 4.77) / 2, + "S36": (0.0129 + 0.0187) / 2, + }, + "Cl": { + "Cl35": (75.5 + 76.1) / 2, + "Cl37": (23.9 + 24.5) / 2, + }, + "Ar": { + "Ar36": 0.3336, + "Ar38": 0.0629, + "Ar40": 99.6035, + }, + "K": { + "K39": 93.2581, + "K40": 0.0117, + "K41": 6.7302, + }, + "Ca": { + "Ca40": 96.941, + "Ca42": 0.647, + "Ca43": 0.135, + "Ca44": 2.086, + "Ca46": 0.004, + "Ca48": 0.187, + }, + "Sc": { + "Sc45": 100.0, + }, + "Ti": { + "Ti46": 8.25, + "Ti47": 7.44, + "Ti48": 73.72, + "Ti49": 5.41, + "Ti50": 5.18, + }, + "V": { + "V50": 0.25, + "V51": 99.75, + }, + "Cr": { + "Cr50": 4.345, + "Cr52": 83.789, + "Cr53": 9.501, + "Cr54": 2.365, + }, + "Mn": { + "Mn55": 100.0, + }, + "Fe": { + "Fe54": 5.845, + "Fe56": 91.754, + "Fe57": 2.119, + "Fe58": 0.282, + }, + "Co": { + "Co59": 100.0, + }, + "Ni": { + "Ni58": 68.0769, + "Ni60": 26.2231, + "Ni61": 1.1399, + "Ni62": 3.6345, + "Ni64": 0.9256, + }, + "Cu": { + "Cu63": 69.15, + "Cu65": 30.85, + }, + "Zn": { + "Zn64": 49.17, + "Zn66": 27.73, + "Zn67": 4.04, + "Zn68": 18.45, + "Zn70": 0.61, + }, + "Ga": { + "Ga69": 60.108, + "Ga71": 39.892, + }, + "Ge": { + "Ge70": 20.52, + "Ge72": 27.45, + "Ge73": 7.76, + "Ge74": 36.52, + "Ge76": 7.75, + }, + "As": { + "As75": 100.0, + }, + "Se": { + "Se74": 0.86, + "Se76": 9.23, + "Se77": 7.6, + "Se78": 23.69, + "Se80": 49.8, + "Se82": 8.82, + }, + "Br": { + "Br79": (50.5 + 50.8) / 2, + "Br81": (49.2 + 49.5) / 2, + }, + "Kr": { + "Kr78": 0.355, + "Kr80": 2.286, + "Kr82": 11.593, + "Kr83": 11.5, + "Kr84": 56.987, + "Kr86": 17.279, + }, + "Rb": { + "Rb85": 72.17, + "Rb87": 27.83, + }, + "Sr": { + "Sr84": 0.56, + "Sr86": 9.86, + "Sr87": 7.0, + "Sr88": 82.58, + }, + "Y": { + "Y89": 100.0, + }, + "Zr": { + "Zr90": 51.45, + "Zr91": 11.22, + "Zr92": 17.15, + "Zr94": 17.38, + "Zr96": 2.8, + }, + "Nb": { + "Nb93": 100.0, + }, + "Mo": { + "Mo92": 14.649, + "Mo94": 9.187, + "Mo95": 15.873, + "Mo96": 16.673, + "Mo97": 9.582, + "Mo98": 24.292, + "Mo100": 9.744, + }, + "Ru": { + "Ru96": 5.54, + "Ru98": 1.87, + "Ru99": 12.76, + "Ru100": 12.6, + "Ru101": 17.06, + "Ru102": 31.55, + "Ru104": 18.62, + }, + "Rh": { + "Rh103": 100.0, + }, + "Pd": { + "Pd102": 1.02, + "Pd104": 11.14, + "Pd105": 22.33, + "Pd106": 27.33, + "Pd108": 26.46, + "Pd110": 11.72, + }, + "Ag": { + "Ag107": 51.839, + "Ag109": 48.161, + }, + "Cd": { + "Cd106": 1.245, + "Cd108": 0.888, + "Cd110": 12.47, + "Cd111": 12.795, + "Cd112": 24.109, + "Cd113": 12.227, + "Cd114": 28.754, + "Cd116": 7.512, + }, + "In": { + "In113": 4.281, + "In115": 95.719, + }, + "Sn": { + "Sn112": 0.97, + "Sn114": 0.66, + "Sn115": 0.34, + "Sn116": 14.54, + "Sn117": 7.68, + "Sn118": 24.22, + "Sn119": 8.59, + "Sn120": 32.58, + "Sn122": 4.63, + "Sn124": 5.79, + }, + "Sb": { + "Sb121": 57.21, + "Sb123": 42.79, + }, + "Te": { + "Te120": 0.09, + "Te122": 2.55, + "Te123": 0.89, + "Te124": 4.74, + "Te125": 7.07, + "Te126": 8.84, + "Te128": 31.74, + "Te130": 34.08, + }, + "I": { + "I127": 100.0, + }, + "Xe": { + "Xe124": 0.095, + "Xe126": 0.089, + "Xe128": 1.91, + "Xe129": 26.4, + "Xe130": 4.071, + "Xe131": 21.232, + "Xe132": 26.909, + "Xe134": 10.436, + "Xe136": 8.857, + }, + "Cs": { + "Cs133": 100.0, + }, + "Ba": { + "Ba130": 0.11, + "Ba132": 0.1, + "Ba134": 2.42, + "Ba135": 6.59, + "Ba136": 7.85, + "Ba137": 11.23, + "Ba138": 71.7, + }, + "La": { + "La138": 0.08881, + "La139": 99.91119, + }, + "Ce": { + "Ce136": 0.186, + "Ce138": 0.251, + "Ce140": 88.449, + "Ce142": 11.114, + }, + "Pr": { + "Pr141": 100.0, + }, + "Nd": { + "Nd142": 27.153, + "Nd143": 12.173, + "Nd144": 23.798, + "Nd145": 8.293, + "Nd146": 17.189, + "Nd148": 5.756, + "Nd150": 5.638, + }, + "Sm": { + "Sm144": 3.08, + "Sm147": 15.0, + "Sm149": 13.82, + "Sm150": 7.37, + "Sm152": 26.74, + "Sm154": 22.74, + }, + "Eu": { + "Eu151": 47.81, + "Eu153": 52.19, + }, + "Gd": { + "Gd152": 0.2, + "Gd154": 2.18, + "Gd155": 14.8, + "Gd156": 20.47, + "Gd157": 15.65, + "Gd158": 24.84, + "Gd160": 21.86, + }, + "Tb": { + "Tb159": 100.0, + }, + "Dy": { + "Dy156": 0.056, + "Dy158": 0.095, + "Dy160": 2.329, + "Dy161": 18.889, + "Dy162": 25.475, + "Dy163": 24.896, + "Dy164": 28.26, + }, + "Ho": { + "Ho165": 100.0, + }, + "Er": { + "Er162": 0.139, + "Er164": 1.601, + "Er166": 33.503, + "Er167": 22.869, + "Er168": 26.978, + "Er170": 14.91, + }, + "Tm": { + "Tm169": 100.0, + }, + "Yb": { + "Yb168": 0.123, + "Yb170": 2.982, + "Yb171": 14.086, + "Yb172": 21.686, + "Yb173": 16.103, + "Yb174": 32.025, + "Yb176": 12.995, + }, + "Lu": { + "Lu175": 97.401, + "Lu176": 2.599, + }, + "Hf": { + "Hf174": 0.16, + "Hf176": 5.26, + "Hf177": 18.6, + "Hf178": 27.28, + "Hf179": 13.62, + "Hf180": 35.08, + }, + "Ta": { + "Ta180": 0.01201, + "Ta181": 99.98799, + }, + "W": { + "W180": 0.12, + "W182": 26.5, + "W183": 14.31, + "W184": 30.64, + "W186": 28.43, + }, + "Re": { + "Re185": 37.4, + "Re187": 62.6, + }, + "Os": { + "Os184": 0.02, + "Os186": 1.59, + "Os187": 1.97, + "Os188": 13.24, + "Os189": 16.15, + "Os190": 26.26, + "Os192": 40.78, + }, + "Ir": { + "Ir191": 37.3, + "Ir193": 62.7, + }, + "Pt": { + "Pt190": 0.012, + "Pt192": 0.782, + "Pt194": 32.864, + "Pt195": 33.77, + "Pt196": 25.21, + "Pt198": 7.356, + }, + "Au": { + "Au197": 100.0, + }, + "Hg": { + "Hg196": 0.15, + "Hg198": 10.04, + "Hg199": 16.94, + "Hg200": 23.14, + "Hg201": 13.17, + "Hg202": 29.74, + "Hg204": 6.82, + }, + "Tl": { + "Tl203": (29.44 + 29.59) / 2, + "Tl205": (70.41 + 70.56) / 2, + }, + "Pb": { + "Pb204": 1.4, + "Pb206": 24.1, + "Pb207": 22.1, + "Pb208": 52.4, + }, + "Bi": { + "Bi209": 100.0, + }, + "Th": { + "Th230": 0.02, + "Th232": 99.98, + }, + "Pa": { + "Pa231": 100.0, + }, + "U": { + "U234": 0.0054, + "U235": 0.7204, + "U238": 99.2742, + }, +} def move_object(object_, velocities, durations): object_.moving = True object_.N_move = len(durations) + 1 diff --git a/mcdc/transport/physics/__init__.py b/mcdc/transport/physics/__init__.py index 73dc04af..66133009 100644 --- a/mcdc/transport/physics/__init__.py +++ b/mcdc/transport/physics/__init__.py @@ -5,4 +5,5 @@ collision_distance, collision, ) +import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron diff --git a/mcdc/transport/physics/electron/__init__.py b/mcdc/transport/physics/electron/__init__.py new file mode 100644 index 00000000..7a0a5f2d --- /dev/null +++ b/mcdc/transport/physics/electron/__init__.py @@ -0,0 +1,6 @@ +from .interface import ( + particle_speed, + macro_xs, + collision, +) +import mcdc.transport.physics.electron.native as native diff --git a/mcdc/transport/physics/electron/interface.py b/mcdc/transport/physics/electron/interface.py new file mode 100644 index 00000000..8fc991ff --- /dev/null +++ b/mcdc/transport/physics/electron/interface.py @@ -0,0 +1,35 @@ +from numba import njit + +#### + +import mcdc.transport.physics.electron.native as native + + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container, mcdc, data): + return native.particle_speed(particle_container) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, mcdc, data): + return native.macro_xs(reaction_type, particle_container, mcdc, data) + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, prog, data): + return native.collision(particle_container, prog, data) diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py new file mode 100644 index 00000000..bbef71b5 --- /dev/null +++ b/mcdc/transport/physics/electron/native.py @@ -0,0 +1,558 @@ +import math +import numpy as np + +from mcdc.mcdc_get import reaction +from numba import njit + +#### + +import mcdc.code_factory.adapt as adapt +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 + +from mcdc.constant import ( + LIGHT_SPEED, + PI, + FINE_STRUCTURE_CONSTANT, + ELECTRON_MASS, + ELECTRON_CUTOFF_ENERGY, + REACTION_TOTAL, + REACTION_ELECTRON_IONIZATION, + REACTION_ELECTRON_ELASTIC_SCATTERING, + REACTION_ELECTRON_BREMSSTRAHLUNG, + REACTION_ELECTRON_EXCITATION, +) +from mcdc.transport.data import evaluate_data +from mcdc.transport.distribution import ( + sample_distribution, + sample_multi_table, +) +from mcdc.transport.physics.util import evaluate_xs_energy_grid, scatter_direction +from mcdc.transport.util import linear_interpolation + + +# ====================================================================================== +# Particle attributes +# ====================================================================================== + + +@njit +def particle_speed(particle_container): + particle = particle_container[0] + E = particle["E"] + mass = ELECTRON_MASS + return LIGHT_SPEED * math.sqrt(E * (E + 2.0 * mass)) / (E + mass) + + +# ====================================================================================== +# Material properties +# ====================================================================================== + + +@njit +def macro_xs(reaction_type, particle_container, mcdc, data): + particle = particle_container[0] + material = mcdc["native_materials"][particle["material_ID"]] + E = particle["E"] + + total = 0.0 + for i in range(material["N_element"]): + element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) + element = mcdc["elements"][element_ID] + + element_density = mcdc_get.native_material.element_densities(i, material, data) + xs = total_micro_xs(reaction_type, E, element, data) + total += element_density * xs + + return total + + +@njit +def total_micro_xs(reaction_type, E, element, data): + idx, E0, E1 = evaluate_xs_energy_grid(E, element, data) + if reaction_type == REACTION_TOTAL: + xs0 = mcdc_get.element.total_xs(idx, element, data) + xs1 = mcdc_get.element.total_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_IONIZATION: + xs0 = mcdc_get.element.ionization_xs(idx, element, data) + xs1 = mcdc_get.element.ionization_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_ELASTIC_SCATTERING: + xs0 = mcdc_get.element.elastic_xs(idx, element, data) + xs1 = mcdc_get.element.elastic_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_EXCITATION: + xs0 = mcdc_get.element.excitation_xs(idx, element, data) + xs1 = mcdc_get.element.excitation_xs(idx + 1, element, data) + elif reaction_type == REACTION_ELECTRON_BREMSSTRAHLUNG: + xs0 = mcdc_get.element.bremsstrahlung_xs(idx, element, data) + xs1 = mcdc_get.element.bremsstrahlung_xs(idx + 1, element, data) + return linear_interpolation(E, E0, E1, xs0, xs1) + + +@njit +def reaction_micro_xs(E, reaction_base, element, data): + idx, E0, E1 = evaluate_xs_energy_grid(E, element, data) + + # Apply offset + offset = reaction_base["xs_offset_"] + if idx < offset: + return 0.0 + idx -= offset + + xs0 = mcdc_get.reaction.xs(idx, reaction_base, data) + xs1 = mcdc_get.reaction.xs(idx + 1, reaction_base, data) + return linear_interpolation(E, E0, E1, xs0, xs1) + + +# ====================================================================================== +# Collision +# ====================================================================================== + + +@njit +def collision(particle_container, prog, data): + mcdc = adapt.mcdc_global(prog) + + particle = particle_container[0] + material = mcdc["native_materials"][particle["material_ID"]] + + # Particle properties + E = particle["E"] + + # Check for cutoff energy + if E <= ELECTRON_CUTOFF_ENERGY: + particle["alive"] = False + edep = E + return edep + + # ================================================================================== + # Sample colliding element + # ================================================================================== + + SigmaT = macro_xs(REACTION_TOTAL, particle_container, mcdc, data) + + # Loop over elements in the material + xi = rng.lcg(particle_container) * SigmaT + total = 0.0 + for i in range(material["N_element"]): + element_ID = int(mcdc_get.native_material.element_IDs(i, material, data)) + element = mcdc["elements"][element_ID] + + element_density = mcdc_get.native_material.element_densities(i, material, data) + sigmaT = total_micro_xs(REACTION_TOTAL, E, element, data) + + total += element_density * sigmaT + + if total > xi: + break + + # ================================================================================== + # Sample and perform reaction + # ================================================================================== + # TODO: The switches for the different reactions have a common code pattern. + # Modularizing the pattern would improve maintainability. + # Note that this also applies to other particle physics. + + sigma_ionization = total_micro_xs(REACTION_ELECTRON_IONIZATION, E, element, data) + sigma_elastic = total_micro_xs( + REACTION_ELECTRON_ELASTIC_SCATTERING, E, element, data + ) + sigma_bremsstrahlung = total_micro_xs( + REACTION_ELECTRON_BREMSSTRAHLUNG, E, element, data + ) + sigma_excitation = total_micro_xs(REACTION_ELECTRON_EXCITATION, E, element, data) + + xi = rng.lcg(particle_container) * sigmaT + total = 0.0 + + # Ionization + total += sigma_ionization + if xi < total: + total -= sigma_ionization + for i in range(element["N_ionization_reaction"]): + reaction_ID = int( + mcdc_get.element.ionization_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_ionization_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return ionization(reaction, particle_container, element, prog, data) + + # Elastic scattering + total += sigma_elastic + if xi < total: + total -= sigma_elastic + for i in range(element["N_elastic_scattering_reaction"]): + reaction_ID = int( + mcdc_get.element.elastic_scattering_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_elastic_scattering_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return elastic_scattering( + reaction, particle_container, element, prog, data + ) + + # Bremsstrahlung + total += sigma_bremsstrahlung + if xi < total: + total -= sigma_bremsstrahlung + for i in range(element["N_bremsstrahlung_reaction"]): + reaction_ID = int( + mcdc_get.element.bremsstrahlung_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_bremsstrahlung_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return bremsstrahlung(reaction, particle_container, element, prog, data) + + # Excitation + total += sigma_excitation + if xi < total: + total -= sigma_excitation + + for i in range(element["N_excitation_reaction"]): + reaction_ID = int( + mcdc_get.element.excitation_reaction_IDs(i, element, data) + ) + reaction = mcdc["electron_excitation_reactions"][reaction_ID] + reaction_base_ID = reaction["parent_ID"] + reaction_base = mcdc["reactions"][reaction_base_ID] + total += reaction_micro_xs(E, reaction_base, element, data) + + if xi < total: + return excitation(reaction, particle_container, element, prog, data) + + +# ====================================================================================== +# Elastic scattering +# ====================================================================================== + + +@njit +def elastic_scattering(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + # ------------------------------------------------------------------------- + # Total elastic xs + # ------------------------------------------------------------------------- + reaction_base_ID = int(reaction["parent_ID"]) + reaction_base = mcdc["reactions"][reaction_base_ID] + xs_total = reaction_micro_xs(E, reaction_base, element, data) + + # If large-angle, xs from data table (EEDL MT=525) + xs_large = 0.0 + xs_large_ID = int(reaction["xs_large_ID"]) + + if xs_large_ID >= 0: + xs_large = elastic_large_xs(E, reaction, mcdc, data) + + # Important to check because of numerical issues + if xs_large < 0.0: + xs_large = 0.0 + if xs_large > xs_total: + xs_large = xs_total + + prob_large = xs_large / xs_total + mu_cut = float(reaction["mu_cut"]) + + xi = rng.lcg(particle_container) + + if xi < prob_large: + # --------------------------------------------------------------------- + # Large-angle elastic scattering + # --------------------------------------------------------------------- + + # Sample mu from EEDL tabulated distribution + mu_table_ID = int(reaction["mu_table_ID"]) + multi_table = mcdc["multi_table_distributions"][mu_table_ID] + mu0 = sample_multi_table(E, particle_container, multi_table, data) + + else: + # --------------------------------------------------------------------- + # Small-angle elastic scattering (Coulomb tail sampling) + # --------------------------------------------------------------------- + + Z = int(element["atomic_number"]) + mu0 = sample_small_angle_mu_coulomb(E, Z, particle_container, mu_cut) + + # Update direction + azi = 2.0 * PI * rng.lcg(particle_container) + # Current direction + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) + + particle["ux"] = ux_new + particle["uy"] = uy_new + particle["uz"] = uz_new + + return 0.0 + + +@njit +def compute_scattering_eta(E, Z): + pc = math.sqrt(E * (E + 2 * ELECTRON_MASS)) + beta = pc / (E + ELECTRON_MASS) + tau = E / ELECTRON_MASS + + r = (FINE_STRUCTURE_CONSTANT * ELECTRON_MASS) / (0.885 * pc) + z_sq = float(Z) ** (2.0 / 3.0) + bracket = 1.13 + 3.76 * ((FINE_STRUCTURE_CONSTANT * float(Z)) / beta) ** 2 + rel = math.sqrt(tau / (tau + 1.0)) + + return 0.25 * (r * r) * z_sq * bracket * rel + + +@njit +def sample_small_angle_mu_coulomb(E, Z, rng_state, mu_cut): + eta = compute_scattering_eta(E, Z) + + x_cut = 1.0 - mu_cut + u = rng.lcg(rng_state) + + denom = (1.0 / eta) - (1.0 / (eta + x_cut)) + inv = (1.0 / eta) - u * denom + x = (1.0 / inv) - eta + + return 1.0 - x + + +@njit +def elastic_large_xs(E, reaction, mcdc, data): + xs_id = int(reaction["xs_large_ID"]) + data_base = mcdc["data"][xs_id] + return evaluate_data(E, data_base, mcdc, data) + + +# ============================================================================== +# Excitation (photon not tracked - all energy deposited) +# ============================================================================== + + +@njit +def excitation(reaction, particle_container, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + dE = evaluate_eloss(E, reaction, mcdc, data) + + # Calculate outgoing energy + E_out = E - dE + + # Check for cutoff + if E_out <= ELECTRON_CUTOFF_ENERGY: + # deposit remaining and kill + edep = E + particle["E"] = 0.0 + particle["alive"] = False + return edep + + # If above cutoff, just deposit dE + particle["E"] = E_out + edep = dE + return edep + + +@njit +def evaluate_eloss(E, reaction, mcdc, data): + data_base = mcdc["data"][int(reaction["eloss_ID"])] + return evaluate_data(E, data_base, mcdc, data) + + +# ============================================================================== +# Bremsstrahlung (photon not tracked - energy deposit = 0 ) +# ============================================================================== + + +@njit +def bremsstrahlung(reaction, particle_container, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + dE = evaluate_eloss(E, reaction, mcdc, data) + + E_out = E - dE + # Check for cutoff + if E_out <= ELECTRON_CUTOFF_ENERGY: + # deposit remaining and kill + edep = E_out + particle["E"] = 0.0 + particle["alive"] = False + return edep + + # If above cutoff, allow photon to escape, just update electron energy + particle["E"] = E_out + edep = 0.0 + return edep + + +# ============================================================================== +# Ionization +# ============================================================================== + + +@njit +def ionization(reaction, particle_container, element, prog, data): + mcdc = adapt.mcdc_global(prog) + + # Particle attributes + particle = particle_container[0] + + # Current energy + E = particle["E"] + + # Sample subshell + N = int(reaction["N_subshell"]) + xs_vals = np.empty(N, dtype=np.float64) + + total = 0.0 + for i in range(N): + xs_sub_id = int(reaction["subshell_xs_IDs"][i]) + xs_sub_table = mcdc["data"][xs_sub_id] + xs_sub_i = evaluate_data(E, xs_sub_table, mcdc, data) + xs_vals[i] = xs_sub_i + total += xs_sub_i + + xi = rng.lcg(particle_container) * total + total_acc = 0.0 + chosen = 0 + for i in range(N): + total_acc += xs_vals[i] + if total_acc >= xi: + chosen = i + break + + # Binding energy (from element, shared across reactions) + B = mcdc_get.element.ionization_subshell_binding_energy(chosen, element, data) + # Deposit all energy if below binding energy + if E <= B: + edep = E + particle["alive"] = False + particle["E"] = 0.0 + return edep + + # Sample secondary energy T_delta from distribution + dist_ID = int(reaction["subshell_product_IDs"][chosen]) + dist_base = mcdc["distributions"][dist_ID] + T_delta = sample_distribution( + E, dist_base, particle_container, mcdc, data, scale=True + ) + + # Primary outgoing energy + E_out = E - B - T_delta + particle["E"] = E_out + + edep = B + + primary_alive_after = True + if E_out <= ELECTRON_CUTOFF_ENERGY: + edep += E_out + particle["E"] = 0.0 + particle["alive"] = False + primary_alive_after = False + + if T_delta <= ELECTRON_CUTOFF_ENERGY: + edep += T_delta + return edep + + # Sample delta direction + ux_delta, uy_delta, uz_delta = sample_delta_direction( + T_delta, E, particle_container + ) + + # Momentum conservation if primary alive + if primary_alive_after: + p_before = math.sqrt(E * (E + 2 * ELECTRON_MASS)) + p_delta = math.sqrt(T_delta * (T_delta + 2 * ELECTRON_MASS)) + + ux_before = particle["ux"] + uy_before = particle["uy"] + uz_before = particle["uz"] + + # Momentum vectors after collision + px_after = p_before * ux_before - p_delta * ux_delta + py_after = p_before * uy_before - p_delta * uy_delta + pz_after = p_before * uz_before - p_delta * uz_delta + + # Normalize and set primary's new direction + norm_sq = px_after * px_after + py_after * py_after + pz_after * pz_after + if norm_sq > 0.0: + norm = math.sqrt(norm_sq) + particle["ux"] = px_after / norm + particle["uy"] = py_after / norm + particle["uz"] = pz_after / norm + + # Add secondary particle to bank + particle_container_new = np.zeros(1, type_.particle_data) + particle_new = particle_container_new[0] + particle_module.copy_as_child(particle_container_new, particle_container) + + particle_new["E"] = T_delta + particle_new["ux"] = ux_delta + particle_new["uy"] = uy_delta + particle_new["uz"] = uz_delta + particle_new["w"] = particle["w"] + + particle_bank_module.add_active(particle_container_new, prog) + return edep + + +@njit +def compute_mu_delta(T_delta, T_prim): + me = ELECTRON_MASS + pd = math.sqrt(T_delta * (T_delta + 2.0 * me)) + pp = math.sqrt(T_prim * (T_prim + 2.0 * me)) + mu = (T_delta * (T_prim + 2.0 * me)) / (pd * pp) + # Check in case of numerical issues + if mu < -1.0: + mu = -1.0 + if mu > 1.0: + mu = 1.0 + + return mu + + +@njit +def sample_delta_direction(T_delta, T_prim, particle_container): + mu = compute_mu_delta(T_delta, T_prim) + phi = 2.0 * PI * rng.lcg(particle_container) + s = math.sqrt(1.0 - mu * mu) + ux = s * math.cos(phi) + uy = s * math.sin(phi) + uz = mu + return ux, uy, uz + + +@njit +def subshell_xs(E, xs_id, mcdc, data): + data_base = mcdc["data"][int(xs_id)] + return evaluate_data(E, data_base, mcdc, data) diff --git a/mcdc/transport/physics/interface.py b/mcdc/transport/physics/interface.py index 34da6f21..a9949427 100644 --- a/mcdc/transport/physics/interface.py +++ b/mcdc/transport/physics/interface.py @@ -5,6 +5,7 @@ ### import mcdc.transport.rng as rng +import mcdc.transport.physics.electron as electron import mcdc.transport.physics.neutron as neutron from mcdc.constant import * @@ -19,6 +20,8 @@ def particle_speed(particle_container, mcdc, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.particle_speed(particle_container, mcdc, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + return electron.particle_speed(particle_container, mcdc, data) return -1.0 @@ -32,17 +35,22 @@ def macro_xs(reaction_type, particle_container, mcdc, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.macro_xs(reaction_type, particle_container, mcdc, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + return electron.macro_xs(reaction_type, particle_container, mcdc, data) return -1.0 @njit def neutron_production_xs(reaction_type, particle_container, mcdc, data): + """ + This is currently only used for neutron k-eigenvalue simulation + """ particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: return neutron.neutron_production_xs( reaction_type, particle_container, mcdc, data ) - return -1.0 + return 0.0 # ====================================================================================== @@ -69,4 +77,6 @@ def collision_distance(particle_container, mcdc, data): def collision(particle_container, mcdc, data): particle = particle_container[0] if particle["particle_type"] == PARTICLE_NEUTRON: - neutron.collision(particle_container, mcdc, data) + neutron.collision(particle_container, prog, data) + elif particle["particle_type"] == PARTICLE_ELECTRON: + electron.collision(particle_container, prog, data) diff --git a/mcdc/transport/physics/neutron/native.py b/mcdc/transport/physics/neutron/native.py index 241e60a6..5f3a918c 100644 --- a/mcdc/transport/physics/neutron/native.py +++ b/mcdc/transport/physics/neutron/native.py @@ -232,6 +232,7 @@ def collision(particle_container, mcdc, data): particle["w"] *= (SigmaT - SigmaC) / SigmaT SigmaT -= SigmaC + # Loop over nuclides in the material xi = rng.lcg(particle_container) * SigmaT total = 0.0 for i in range(material["N_nuclide"]): @@ -246,8 +247,7 @@ def collision(particle_container, mcdc, data): particle["w"] *= (sigmaT - sigmaC) / sigmaT sigmaT -= sigmaC - SigmaT_nuclide = nuclide_density * sigmaT - total += SigmaT_nuclide + total += nuclide_density * sigmaT if total > xi: break @@ -255,6 +255,9 @@ def collision(particle_container, mcdc, data): # ================================================================================== # Sample and perform reaction # ================================================================================== + # TODO: The switches for the different reactions have a common code pattern. + # Modularizing the pattern would improve maintainability. + # Note that this also applies to other particle physics. sigma_elastic = total_micro_xs( REACTION_NEUTRON_ELASTIC_SCATTERING, E, nuclide, data @@ -294,7 +297,6 @@ def collision(particle_container, mcdc, data): total += sigma_inelastic if xi < total: total -= sigma_inelastic - for i in range(nuclide["N_inelastic_scattering_reaction"]): reaction_ID = int( mcdc_get.nuclide.inelastic_scattering_reaction_IDs(i, nuclide, data) @@ -302,8 +304,7 @@ def collision(particle_container, mcdc, data): reaction = mcdc["neutron_inelastic_scattering_reactions"][reaction_ID] reaction_base_ID = reaction["parent_ID"] reaction_base = mcdc["reactions"][reaction_base_ID] - xs = reaction_micro_xs(E, reaction_base, nuclide, data) - total += xs + total += reaction_micro_xs(E, reaction_base, nuclide, data) if xi < total: inelastic_scattering(reaction, particle_container, nuclide, mcdc, data) return diff --git a/test/regression/electron_tests/lockwood/experimental_conf_Al.py b/test/regression/electron_tests/lockwood/experimental_conf_Al.py new file mode 100644 index 00000000..00663d6d --- /dev/null +++ b/test/regression/electron_tests/lockwood/experimental_conf_Al.py @@ -0,0 +1,133 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +# Set the XS library directory +os.environ["MCDC_XSLIB"] = os.path.join(os.path.dirname(os.path.dirname(os.getcwd())), "dummy_data") + +# ============================================================================= +# Set problem parameters +# ============================================================================= +MATERIAL_SYMBOL = "Al" +CSDA_RANGE = 0.569 # g/cm2 +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Other parameters (don't change) +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +TINY = 1e-12 +VOID_SCALE = 1.0e-30 +VOID_DENSITY = MAT_DENSITY_ATOMS_PER_BARN_CM * VOID_SCALE + +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) + +# ============================================================================= +# Exactly the same computational configuration as MCNP for Lockwood exp (from the MCNP VV slides) +# https://mcnp.lanl.gov/pdf_files/TechReport_2024_LANL_LA-UR-23-32743Rev.1_Kulesza.pdf +# ============================================================================= +# pz 0.0 +# pz 0.0233 - front foil thickness = 0.0233 cm +# pz 0.1233 - +0.1 cm void gap +# pz 0.12517037037 - calorimeter foil thickness = 0.00187037037 cm (5.05e-3/2.7) +# pz 0.22517037037 - +0.1 cm void gap +# pz 5.0 - thick Al plate end +# +# cz 100.0 - radius 100 cm boundary +z0 = 0.0 +z1 = 0.0233 +z2 = 0.1233 +z3 = 0.12517037037037037 +z4 = 0.22517037037037037 +z5 = 5.0 + +R_CYL = 100.0 + +# ============================================================================= +# Set materials +# ============================================================================= +mat_al = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +mat_void = mcdc.Material( + element_composition={MATERIAL_SYMBOL: VOID_DENSITY} +) + +# ============================================================================= +# Surfaces +# ============================================================================= +# Z-direction surfaces for layers + +s0 = mcdc.surface("plane-z", z=z0, bc="vacuum") +s1 = mcdc.surface("plane-z", z=z1) +s2 = mcdc.surface("plane-z", z=z2) +s3 = mcdc.surface("plane-z", z=z3) +s4 = mcdc.surface("plane-z", z=z4) +s5 = mcdc.surface("plane-z", z=z5, bc="vacuum") + +# ============================================================================= +# Cells +# ============================================================================= + +# Front foil (Al): z0 < z < z1 +mcdc.cell(+s0 & -s1, mat_al) + +# Void gap (almost vacuum): z1 < z < z2 +mcdc.cell(+s1 & -s2, mat_void) + +# Calorimeter foil (Al): z2 < z < z3 +mcdc.cell(+s2 & -s3, mat_al) + +# Void gap: z3 < z < z4 +mcdc.cell(+s3 & -s4, mat_void) + +# Thick Al plate: z4 < z < z5 +mcdc.cell(+s4 & -s5, mat_al) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0, 0 degree incidence +theta = 0.0 +mcdc.source( + z=[z0 + 1e-8, z0 + 1e-8], + particle_type="electron", + energy=np.array([[1e6 - 1.0, 1e6 + 1.0], [0.5, 0.5]]), + direction=[math.sin(theta), TINY, math.cos(theta)], +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(z4, z4+L, 101) + +mcdc.tally.mesh_tally( + scores=["edep", "flux"], + z=z_bins +) + +mcdc.tally.surface_tally(s0, scores=["net-current"]) +mcdc.tally.surface_tally(s5, scores=["net-current"]) + + +# ============================================================================= +# Settings and run +# ============================================================================= +settings = mcdc.Settings( + N_particle=100, + active_bank_buffer=10000 +) + +settings.save_input_deck = True +settings.output_name = f"lockwood_exp_setup_{datetime.now():%Y%m%d_%H%M%S}" +settings.debug_energy = False +settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/test/regression/electron_tests/lockwood/slab_Al.py b/test/regression/electron_tests/lockwood/slab_Al.py new file mode 100644 index 00000000..607af8dc --- /dev/null +++ b/test/regression/electron_tests/lockwood/slab_Al.py @@ -0,0 +1,101 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +# Set the XS library directory +os.environ["MCDC_XSLIB"] = os.path.join(os.path.dirname(os.path.dirname(os.getcwd())), "dummy_data") + +# ============================================================================= +# Set problem parameters +# ============================================================================= +MATERIAL_SYMBOL = "Al" +CSDA_RANGE = 0.569 # g/cm2 +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Other parameters (don't change) +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +TINY = 1e-12 +VOID_SCALE = 1.0e-30 +VOID_DENSITY = MAT_DENSITY_ATOMS_PER_BARN_CM * VOID_SCALE +M = 20 + +L = CSDA_RANGE / RHO_G_CM3 * M # cm +N_LAYERS = int(L / dz) + +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.surface("plane-z", z=0, bc="vacuum") +s2 = mcdc.surface("plane-z", z=L, bc="vacuum") + +mcdc.cell(+s1 & -s2, mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 +theta = 0 +z0 = 0.0 + +mcdc.source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[1e6 - 1, 1e6 + 1], [0.5, 0.5]]), + direction=[math.sin(theta), 0.0+TINY, math.cos(theta)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L/M, N_LAYERS + 1) + +mcdc.tally.mesh_tally( + scores=["edep", "flux"], + z=z_bins +) + +mcdc.tally.surface_tally(s1, scores=["net-current"]) +mcdc.tally.surface_tally(s2, scores=["net-current"]) + + +# ============================================================================= +# Settings and run +# ============================================================================= +settings = mcdc.Settings( + N_particle=100, + active_bank_buffer=10000 +) + +settings.save_input_deck = True +settings.output_name = f"lockwood_output_{datetime.now():%Y%m%d_%H%M%S}" +settings.debug_energy = False +settings.use_progress_bar = True + + +mcdc.run() \ No newline at end of file diff --git a/test/regression/electron_tests/references.md b/test/regression/electron_tests/references.md new file mode 100644 index 00000000..22e07924 --- /dev/null +++ b/test/regression/electron_tests/references.md @@ -0,0 +1,7 @@ +# References used during this work. + +1 - [Calorimetric Measurement of Electron Energy Deposition In Extended Media Theory vs Experiment - Lockwood](https://www.osti.gov/servlets/purl/5516784) +1 - [Validation of MCNP6 for Electron Energy Deposition in Extended Media](https://mcnp.lanl.gov/pdf_files/TechReport_2015_LANL_LA-UR-15-25134_DixonHughes.pdf) +2 - [An Electron/Photon/Relaxation Data Library for MCNP6](https://nucleardata.lanl.gov/files/la-ur-13-27377.pdf) +3 - [Electron/Photon/Relaxation Data Library—eprdata14](https://nucleardata.lanl.gov/ace/eprdata14/) +4 - [MCNP6.3 Electron Energy Deposition Validation with the Lockwood Experiments](https://mcnp.lanl.gov/pdf_files/TechReport_2024_LANL_LA-UR-23-32743Rev.1_Kulesza.pdf) \ No newline at end of file