diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index 3074001..374a9a5 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -28,7 +28,6 @@ jobs: - name: install package run: | conda install -c conda-forge openmc - # conda install -c cadquery -c conda-forge cadquery=2.1 python setup.py install - name: install packages for tests @@ -45,6 +44,7 @@ jobs: run: | cd examples python processing_2d_mesh_effective_dose_tally.py + python processing_cell_damage-energy_tally.py python processing_cell_effective_dose_tally.py python processing_cell_flux_tally.py python processing_cell_heating_tally.py diff --git a/examples/processing_2d_mesh_effective_dose_tally.py b/examples/processing_2d_mesh_effective_dose_tally.py index 1681d8c..b52e420 100644 --- a/examples/processing_2d_mesh_effective_dose_tally.py +++ b/examples/processing_2d_mesh_effective_dose_tally.py @@ -17,8 +17,7 @@ # scaled from picosievert to sievert result = opp.process_dose_tally( - tally=my_tally, - required_units="sievert cm **2 / simulated_particle" + tally=my_tally, required_units="sievert cm **2 / simulated_particle" ) # the tally result with required units print(result) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py new file mode 100644 index 0000000..e0fede1 --- /dev/null +++ b/examples/processing_cell_damage-energy_tally.py @@ -0,0 +1,68 @@ +import openmc_post_processor as opp +import openmc + +# loads in the statepoint file containing tallies +statepoint = openmc.StatePoint(filepath="statepoint.2.h5") + +# gets one tally from the available tallies +my_tally = statepoint.get_tally(name="2_damage-energy") + + +# returns the tally with base units +result = opp.process_damage_energy_tally( + tally=my_tally, required_units="eV / simulated_particle" +) +print(f"damage-energy base units = {result}", end="\n\n") + + +# returns the tally with scalled based units (MeV instead of eV) +result, std_dev = opp.process_damage_energy_tally( + tally=my_tally, + required_units="MeV / simulated_particle" + # required_units="displacements per atoms" +) +print(f"damage-energy scaled base units = {result}", end="\n\n") + +# # returns the tally with scalled based units (MeV instead of eV) +result, std_dev = opp.process_damage_energy_tally( + tally=my_tally, + required_units="MeV", + source_strength=1, +) +print(f"damage-energy scaled base units = {result}", end="\n\n") + + +my_mat = openmc.Material() +my_mat.add_element("Fe", 1) +my_mat.set_density("g/cm3", 1) + +result, std_dev = opp.process_damage_energy_tally( + tally=my_tally, + required_units="MeV per atom", + source_strength=1, + volume=5, + material=my_mat, +) +print(f"damage-energy scaled base units = {result}", end="\n\n") + +result, std_dev = opp.process_damage_energy_tally( + tally=my_tally, + required_units="displacements", + source_strength=1, + volume=5, + material=my_mat, + energy_per_displacement=10, + recombination_fraction=0.1, +) +print(f"damage-energy = {result}", end="\n\n") + +result, std_dev = opp.process_damage_energy_tally( + tally=my_tally, + required_units="displacements per pulse", + source_strength=0.5, + volume=5, + material=my_mat, + energy_per_displacement=10, + recombination_fraction=0.1, +) +print(f"damage-energy = {result}", end="\n\n") diff --git a/examples/processing_cell_flux_tally.py b/examples/processing_cell_flux_tally.py index d11ed73..2b95d21 100644 --- a/examples/processing_cell_flux_tally.py +++ b/examples/processing_cell_flux_tally.py @@ -11,8 +11,7 @@ # returns the tally with base units result = opp.process_tally( - tally=my_tally, - required_units="centimeter / simulated_particle" + tally=my_tally, required_units="centimeter / simulated_particle" ) print(f"flux base units = {result[0].units}", end="\n\n") diff --git a/examples/processing_cell_heating_tally.py b/examples/processing_cell_heating_tally.py index 1fe5c96..72df1ac 100644 --- a/examples/processing_cell_heating_tally.py +++ b/examples/processing_cell_heating_tally.py @@ -9,10 +9,7 @@ # returns the tally with base units -result = opp.process_tally( - tally=my_tally, - required_units='eV / simulated_particle' -) +result = opp.process_tally(tally=my_tally, required_units="eV / simulated_particle") print(f"heating base units = {result}", end="\n\n") diff --git a/examples/processing_cell_spectra_tally.py b/examples/processing_cell_spectra_tally.py index 9a5d10a..adebbe2 100644 --- a/examples/processing_cell_spectra_tally.py +++ b/examples/processing_cell_spectra_tally.py @@ -29,7 +29,7 @@ tally=my_tally, required_units="centimeter / pulse", required_energy_units="eV", - source_strength=1.3e6 + source_strength=1.3e6, ) print(f"spectra per pulse = {result}", end="\n\n") @@ -40,6 +40,6 @@ tally=my_tally, required_units="centimeter / second", required_energy_units="MeV", - source_strength=1e9 + source_strength=1e9, ) print(f"spectra per pulse = {result}", end="\n\n") diff --git a/examples/processing_mutliple_cell_spectra_tally.py b/examples/processing_mutliple_cell_spectra_tally.py index 9df7da6..f5ee7a4 100644 --- a/examples/processing_mutliple_cell_spectra_tally.py +++ b/examples/processing_mutliple_cell_spectra_tally.py @@ -19,7 +19,7 @@ y_scale="log", filename="combine_spectra_plot_1.html", required_energy_units="eV", - plotting_package='plotly', + plotting_package="plotly", required_units="centimeters / simulated_particle", ) @@ -33,7 +33,7 @@ y_scale="log", filename="combine_spectra_plot_2.html", required_energy_units="MeV", - plotting_package='plotly', + plotting_package="plotly", required_units="centimeter / second", source_strength=1.3e6, ) @@ -47,8 +47,8 @@ y_scale="log", filename="combine_spectra_plot_3.html", required_energy_units="MeV", - plotting_package='plotly', + plotting_package="plotly", required_units="neutrons / second * cm ** -2", source_strength=1.3e7, - volume=100 + volume=100, ) diff --git a/openmc_post_processor/__init__.py b/openmc_post_processor/__init__.py index 1887241..002bf1b 100644 --- a/openmc_post_processor/__init__.py +++ b/openmc_post_processor/__init__.py @@ -7,5 +7,6 @@ process_tally, process_dose_tally, process_spectra_tally, - scale_tally + process_damage_energy_tally, + scale_tally, ) diff --git a/openmc_post_processor/neutronics_units.txt b/openmc_post_processor/neutronics_units.txt index 4e771af..8388bde 100644 --- a/openmc_post_processor/neutronics_units.txt +++ b/openmc_post_processor/neutronics_units.txt @@ -5,4 +5,6 @@ particle = [] neutron = [] photon = [] simulated_particle = [] +atom = [atom] +displacements = [displacements] pulse = [pulse] diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 735dac3..49bfa6c 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -9,10 +9,109 @@ ureg.load_definitions(str(Path(__file__).parent / "neutronics_units.txt")) +def process_damage_energy_tally( + tally, + required_units: str = "eV / simulated_particle", + source_strength: float = None, + volume: float = None, + energy_per_displacement: float = None, + recombination_fraction: float = 0, + material: float = None, +): + """Processes a damage-energy tally converting the tally with default units + obtained during simulation into the user specified units. Can be processed + to obtain damage per atom (DPA) + + Args: + tally: The openmc.Tally object which should be a spectra tally. With a + score of flux or current and an EnergyFunctionFilter + required_units: The units to convert the energy and tally into + source_strength: In some cases the source_strength will be required + to convert the base units into the required units. This optional + argument allows the user to specify the source_strength when needed + volume: In some cases the volume will be required to convert the base + units into the required units. In the case of a regular mesh the + volume is automatically found. This optional argument allows the + user to specify the volume when needed or overwrite the + automatically calculated volume. When finding DPA volume is needed + along with the material to find number of atoms. + energy_per_displacement: the energy required to displace an atom. The + total damage-energy depositied is divided by this value to get + number of atoms displaced. Assumed units are eV + + Returns: + The dpa tally result in the required units + """ + + if check_for_energy_function_filter(tally): + raise ValueError("EnergyFunctionFilter found in a damage-energy tally") + + # checks for user provided base units + base_units = get_tally_units(tally) + + data_frame = tally.get_pandas_dataframe() + + print(f"tally {tally.name} base units {base_units}") + + tally_result = np.array(data_frame["mean"]) + + if recombination_fraction: + if recombination_fraction < 0: + raise ValueError( + f"recombination_fraction can't be smaller than 1. recombination_fraction is {recombination_fraction}" + ) + if recombination_fraction > 1: + raise ValueError( + f"recombination_fraction can't be larger than 1. recombination_fraction is {recombination_fraction}" + ) + + tally_result = tally_result * (1.0 - recombination_fraction) + + tally_result = tally_result * base_units + + if material: + atomic_mass_in_g = material.average_molar_mass * 1.66054e-24 + density_in_g_per_cm3 = material.get_mass_density() + number_of_atoms_per_cm3 = density_in_g_per_cm3 / atomic_mass_in_g + else: + number_of_atoms_per_cm3 = None + print("number_of_atoms_per_cm3", number_of_atoms_per_cm3) + + scaled_tally_result = scale_tally( + tally, + tally_result, + ureg[required_units], + source_strength, + volume, + number_of_atoms_per_cm3, + energy_per_displacement, + ) + + tally_in_required_units = scaled_tally_result.to(required_units) + + if "std. dev." in data_frame.columns.to_list(): + tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units + scaled_tally_std_dev = scale_tally( + tally, + tally_std_dev_base, + ureg[required_units], + source_strength, + volume, + number_of_atoms_per_cm3, + energy_per_displacement, + ) + # if recombination_fraction: + # scaled_tally_std_dev = scaled_tally_std_dev * recombination_fraction + tally_std_dev_in_required_units = scaled_tally_std_dev.to(required_units) + return tally_in_required_units, tally_std_dev_in_required_units + else: + return tally_in_required_units + + def process_spectra_tally( tally, required_units: str = "centimeters / simulated_particle", - required_energy_units: str = 'eV', + required_energy_units: str = "eV", source_strength: float = None, volume: float = None, ) -> tuple: @@ -38,10 +137,10 @@ def process_spectra_tally( """ if not check_for_energy_filter(tally): - raise ValueError('EnergyFilter was not found in spectra tally') + raise ValueError("EnergyFilter was not found in spectra tally") if check_for_energy_function_filter(tally): - raise ValueError('EnergyFunctionFilter was found in spectra tally') + raise ValueError("EnergyFunctionFilter was found in spectra tally") data_frame = tally.get_pandas_dataframe() @@ -58,7 +157,6 @@ def process_spectra_tally( scaled_tally_result = scale_tally( tally, tally_base, - base_units, ureg[required_units], source_strength, volume, @@ -70,14 +168,17 @@ def process_spectra_tally( scaled_tally_std_dev = scale_tally( tally, tally_std_dev_base, - base_units, ureg[required_units], source_strength, volume, ) tally_std_dev_in_required_units = scaled_tally_std_dev.to(required_units) - - return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units + + return ( + energy_in_required_units, + tally_in_required_units, + tally_std_dev_in_required_units, + ) else: @@ -86,7 +187,7 @@ def process_spectra_tally( def process_dose_tally( tally, - required_units: str = 'picosievert cm **2 / simulated_particle', + required_units: str = "picosievert cm **2 / simulated_particle", source_strength: float = None, volume: float = None, ): @@ -111,12 +212,17 @@ def process_dose_tally( """ if not check_for_energy_function_filter(tally): - raise ValueError('EnergyFunctionFilter was not found in dose tally') + raise ValueError("EnergyFunctionFilter was not found in dose tally") # checks for user provided base units base_units = get_tally_units(tally) base_units = ureg.picosievert * ureg.centimeter * base_units + # dose coefficients have pico sievert cm **2 + # flux has cm2 / simulated_particle units + # dose on a surface uses a current score (units of per simulated_particle) and is therefore * area to get pSv / source particle + # dose on a volume uses a flux score (units of cm2 per simulated particle) and therefore gives pSv cm**4 / simulated particle + data_frame = tally.get_pandas_dataframe() print(f"tally {tally.name} base units {base_units}") @@ -126,7 +232,6 @@ def process_dose_tally( scaled_tally_result = scale_tally( tally, tally_result, - base_units, ureg[required_units], source_strength, volume, @@ -138,7 +243,6 @@ def process_dose_tally( scaled_tally_std_dev = scale_tally( tally, tally_std_dev_base, - base_units, ureg[required_units], source_strength, volume, @@ -198,7 +302,6 @@ def process_tally( scaled_tally_result = scale_tally( tally, tally_result, - base_units, ureg[required_units], source_strength, volume, @@ -210,7 +313,6 @@ def process_tally( scaled_tally_std_dev = scale_tally( tally, tally_std_dev_base, - base_units, ureg[required_units], source_strength, volume, @@ -227,13 +329,63 @@ def process_tally( def scale_tally( tally, tally_result, - base_units, - required_units, + required_units: str, source_strength: float, - volume: float + volume: float, + atoms: float = None, + energy_per_displacement: float = None, ): + + # energy_per_displacement time_diff = check_for_dimentionality_difference( - base_units, required_units, "[time]" + tally_result.units, required_units, "[time]" + ) + length_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[length]" + ) + mass_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[mass]" + ) + displacement_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[displacements]" + ) + + if ( + time_diff == -2 + and mass_diff == 1 + and length_diff == 2 + and displacement_diff == -1 + ): + + print("energy per displacement_diff scaling needed (eV)") + if energy_per_displacement: + energy_per_displacement = ( + energy_per_displacement * ureg.electron_volt / ureg["displacements"] + ) + tally_result = tally_result / energy_per_displacement + else: + raise ValueError( + f"energy_per_displacement is required but currently set to {energy_per_displacement}" + ) + + # per_displacement + displacement_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[displacements]" + ) + if displacement_diff == -1: + print("energy per displacement_diff scaling needed (eV)") + if energy_per_displacement: + energy_per_displacement = ( + energy_per_displacement * ureg.electron_volt / ureg["displacements"] + ) + tally_result = tally_result / energy_per_displacement + else: + raise ValueError( + f"energy_per_displacement is required but currently set to {energy_per_displacement}" + ) + + time_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[time]" ) if time_diff != 0: print("time scaling needed (seconds)") @@ -249,7 +401,7 @@ def scale_tally( ) time_diff = check_for_dimentionality_difference( - base_units, required_units, "[pulse]" + tally_result.units, required_units, "[pulse]" ) if time_diff != 0: print("time scaling needed (pulse)") @@ -257,7 +409,7 @@ def scale_tally( source_strength = source_strength * ureg["1 / pulse"] if time_diff == -1: tally_result = tally_result / source_strength - if time_diff == 1: + elif time_diff == 1: tally_result = tally_result * source_strength else: raise ValueError( @@ -265,29 +417,63 @@ def scale_tally( ) length_diff = check_for_dimentionality_difference( - base_units, required_units, "[length]" + tally_result.units, required_units, "[length]" ) if length_diff != 0: print("length scaling needed") if volume: - volume = volume * ureg["centimeter ** 3"] + volume_with_units = volume * ureg["centimeter ** 3"] else: # volume required but not provided so it is found from the mesh - volume = compute_volume_of_voxels(tally) * ureg["centimeter ** 3"] + volume_from_mesh = compute_volume_of_voxels(tally) + + if volume_from_mesh: + volume_with_units = volume_from_mesh * ureg["centimeter ** 3"] + else: + msg = ( + f"A length dimentionality difference of {length_diff} " + f"was detected. However volume is set to {volume} and " + "volume could not be calculated from the mesh. Please " + "specify the volume argument" + ) + raise ValueError(msg) if length_diff == 3: print("dividing by volume") - tally_result = tally_result / volume - if length_diff == -3: + tally_result = tally_result / volume_with_units + elif length_diff == -3: print("multiplying by volume") - tally_result = tally_result * volume + tally_result = tally_result * volume_with_units + atom_diff = check_for_dimentionality_difference( + tally_result.units, required_units, "[atom]" + ) + if atom_diff != 0: + print("atom scaling needed") + if atoms: + atoms = atoms * ureg["atom"] + + if atom_diff == 1: + print("dividing by atom") + tally_result = tally_result / atoms + elif atom_diff == -1: + print("multiplying by atom") + tally_result = tally_result * atoms + + else: + msg = ( + f"atoms is required but currently set to {atoms}. Atoms " + "can be calculated automatically from material and volume " + "inputs" + ) + raise ValueError(msg) return tally_result def compute_volume_of_voxels(tally): - tally_filter = tally.find_filter(filter_type=openmc.MeshFilter) - if tally_filter: + if tally.contains_filter(openmc.MeshFilter): + tally_filter = tally.find_filter(filter_type=openmc.MeshFilter) + mesh = tally_filter.mesh x = abs(mesh.lower_left[0] - mesh.upper_right[0]) / mesh.dimension[0] y = abs(mesh.lower_left[1] - mesh.upper_right[1]) / mesh.dimension[1] @@ -295,7 +481,8 @@ def compute_volume_of_voxels(tally): volume = x * y * z return volume else: - raise ValueError(f"volume could not be obtained from tally {tally}") + print(f"volume of mesh element could not be obtained from tally {tally}") + return False def find_fusion_energy_per_reaction(reactants: str) -> float: @@ -361,29 +548,6 @@ def get_cell_ids_from_tally_filters(tally): return cell_ids -# def get_tally_units_dose(tally): - -# # An EnergyFunctionFilter is exspeted in dose tallies -# units = get_tally_units(tally) - -# units = [ureg.picosievert * ureg.centimeter * units[0]] - -# # check it is a dose tally by looking for a openmc.filter.EnergyFunctionFilter -# for filter in tally.filters: -# if isinstance(filter, openmc.filter.EnergyFunctionFilter): -# print("filter is EnergyFunctionFilter") -# # effective_dose -# # dose coefficients have pico sievert cm **2 -# # flux has cm2 / simulated_particle units -# # dose on a surface uses a current score (units of per simulated_particle) and is therefore * area to get pSv / source particle -# # dose on a volume uses a flux score (units of cm2 per simulated particle) and therefore gives pSv cm**4 / simulated particle - -# return units - -# raise ValueError( -# "units for dose tally can't be found, an EnergyFunctionFilter was not present" -# ) - def check_for_energy_filter(tally): # check it is a spectra tally by looking for a openmc.filter.EnergyFilter @@ -393,10 +557,6 @@ def check_for_energy_filter(tally): return True return False - # raise ValueError( - # "units for spectra tally can't be found, an EnergyFilter was not present" - # ) - def check_for_energy_function_filter(tally): # check for EnergyFunctionFilter which modify the units of the tally @@ -424,16 +584,22 @@ def get_tally_units(tally): # heating units are eV / simulated_particle units = ureg.electron_volt / ureg.simulated_particle + elif tally.scores == ["damage-energy"]: + # damage-energy units are eV / simulated_particle + units = ureg.electron_volt / ureg.simulated_particle + else: - # return [1 / ureg.simulated_particle] - raise ValueError( - "units for tally can't be found, supported tallies are currently limited" + msg = ( + "units for tally can't be found. Tallies that are supported " + "by get_tally_units function are those with scores of current, " + "flux, heating, damage-energy" ) + raise ValueError(msg) return units def check_for_dimentionality_difference(units_1, units_2, unit_to_compare): - units_1_time_power = units_1.dimensionality.get(unit_to_compare) - units_2_time_power = units_2.dimensionality.get(unit_to_compare) - return units_1_time_power - units_2_time_power + units_1_dimentions = units_1.dimensionality.get(unit_to_compare) + units_2_dimentions = units_2.dimensionality.get(unit_to_compare) + return units_1_dimentions - units_2_dimentions diff --git a/setup.py b/setup.py index 8f098fb..8abe81f 100644 --- a/setup.py +++ b/setup.py @@ -14,17 +14,17 @@ url="https://github.com/fusion-energy/openmc_post_processor", packages=setuptools.find_packages(), classifiers=[ - 'Natural Language :: English', - 'Topic :: Scientific/Engineering', - 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.6', - 'Programming Language :: Python :: 3.7', - 'Programming Language :: Python :: 3.8', - 'Programming Language :: Python :: 3.9', - 'License :: OSI Approved :: MIT License', - 'Operating System :: OS Independent', + "Natural Language :: English", + "Topic :: Scientific/Engineering", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.6", + "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", ], - python_requires='>=3.6', + python_requires=">=3.6", package_data={ "openmc_post_processor": [ # "requirements.txt", @@ -33,6 +33,5 @@ "neutronics_units.txt", ] }, - install_requires=[ - "pint" - ]) + install_requires=["pint"], +) diff --git a/tests/create_statepoint_file_for_testing.py b/tests/create_statepoint_file_for_testing.py index 604b463..d28b2e2 100644 --- a/tests/create_statepoint_file_for_testing.py +++ b/tests/create_statepoint_file_for_testing.py @@ -119,14 +119,19 @@ target=2, ) -tally14 = odw.MeshTally2D( +tally14 = odw.CellTally( + tally_type="damage-energy", + target=2, +) + +tally15 = odw.MeshTally2D( tally_type="neutron_effective_dose", plane="xy", mesh_resolution=(10, 5), bounding_box=[(-100, -100, 0), (100, 100, 1)], ) -tally15 = odw.MeshTally3D( +tally16 = odw.MeshTally3D( mesh_resolution=(100, 100, 100), bounding_box=[(-100, -100, 0), (100, 100, 1)], tally_type="neutron_effective_dose", @@ -149,9 +154,11 @@ tally13, tally14, tally15, + tally16, ] ) + settings = odw.FusionSettings() settings.batches = 2 settings.particles = 1000000 diff --git a/tests/create_statepoint_file_for_testing_one_batch.py b/tests/create_statepoint_file_for_testing_one_batch.py index 8dab690..9b33533 100644 --- a/tests/create_statepoint_file_for_testing_one_batch.py +++ b/tests/create_statepoint_file_for_testing_one_batch.py @@ -119,14 +119,19 @@ target=2, ) -tally14 = odw.MeshTally2D( +tally14 = odw.CellTally( + tally_type="damage-energy", + target=2, +) + +tally15 = odw.MeshTally2D( tally_type="neutron_effective_dose", plane="xy", mesh_resolution=(10, 5), bounding_box=[(-100, -100, 0), (100, 100, 1)], ) -tally15 = odw.MeshTally3D( +tally16 = odw.MeshTally3D( mesh_resolution=(100, 100, 100), bounding_box=[(-100, -100, 0), (100, 100, 1)], tally_type="neutron_effective_dose", @@ -149,6 +154,7 @@ tally13, tally14, tally15, + tally16, ] ) diff --git a/tests/create_statepoint_file_for_testing_using_dagmc.py b/tests/create_statepoint_file_for_testing_using_dagmc.py index 13bb8fa..39528ee 100644 --- a/tests/create_statepoint_file_for_testing_using_dagmc.py +++ b/tests/create_statepoint_file_for_testing_using_dagmc.py @@ -101,14 +101,19 @@ target=1, ) -tally14 = odw.MeshTally2D( +tally14 = odw.CellTally( + tally_type="damage-energy", + target=2, +) + +tally15 = odw.MeshTally2D( tally_type="neutron_effective_dose", plane="xy", mesh_resolution=(10, 5), bounding_box=[(-100, -100, 0), (100, 100, 1)], ) -tally15 = odw.MeshTally3D( +tally16 = odw.MeshTally3D( mesh_resolution=(100, 100, 100), bounding_box=[(-100, -100, 0), (100, 100, 1)], tally_type="neutron_effective_dose", @@ -131,9 +136,11 @@ tally13, tally14, tally15, + tally16, ] ) + settings = odw.FusionSettings() settings.batches = 3 settings.particles = 1000000 diff --git a/tests/test_cell_damage_energy.py b/tests/test_cell_damage_energy.py new file mode 100644 index 0000000..b734bd8 --- /dev/null +++ b/tests/test_cell_damage_energy.py @@ -0,0 +1,257 @@ +import unittest + +import openmc_post_processor as opp +import pytest +import openmc + + +class TestUsage(unittest.TestCase): + def setUp(self): + + # loads in the statepoint file containing tallies + statepoint = openmc.StatePoint(filepath="statepoint.2.h5") + self.my_tally = statepoint.get_tally(name="2_damage-energy") + + def test_energy_no_processing(self): + # returns the tally with base units + result = opp.process_damage_energy_tally( + tally=self.my_tally, required_units="eV / simulated_particle" + ) + + assert len(result) == 2 + assert result[0].units == "electron_volt / simulated_particle" + assert result[1].units == "electron_volt / simulated_particle" + assert isinstance(result[0][0].magnitude, float) + assert isinstance(result[1][0].magnitude, float) + + def test_energy_fusion_power_processing(self): + + # returns the tally with scalled based units (MeV instead of eV) + result = opp.process_damage_energy_tally( + source_strength=4.6e17, # neutrons per 1.3MJ pulse + tally=self.my_tally, + required_units="eV / second", + ) + + assert len(result) == 2 + assert result[0].units == "electron_volt / second" + assert result[1].units == "electron_volt / second" + assert isinstance(result[0][0].magnitude, float) + assert isinstance(result[1][0].magnitude, float) + + def test_energy_pulse_processing(self): + + # returns the tally with scalled based units (MeV instead of eV) + result = opp.process_damage_energy_tally( + source_strength=4.6e17, # neutrons per 1.3MJ pulse + tally=self.my_tally, + required_units="eV / pulse", + ) + + assert len(result) == 2 + assert result[0].units == "electron_volt / pulse" + assert result[1].units == "electron_volt / pulse" + + def test_energy_pulse_processing_and_scaling(self): + + # returns the tally with scalled based units (MeV instead of eV) + result = opp.process_damage_energy_tally( + source_strength=4.6e17, # neutrons per 1.3MJ pulse + tally=self.my_tally, + required_units="MeV / pulse", + ) + + assert len(result) == 2 + assert result[0].units == "megaelectron_volt / pulse" + assert result[1].units == "megaelectron_volt / pulse" + + def test_energy_fusion_power_processing_and_scaling(self): + + # returns the tally with scalled based units (MeV instead of eV) + result = opp.process_damage_energy_tally( + source_strength=4.6e17, # neutrons per 1.3MJ pulse + tally=self.my_tally, + required_units="MeV / second", + ) + + assert len(result) == 2 + assert result[0].units == "megaelectron_volt / second" + assert result[1].units == "megaelectron_volt / second" + + def test_energy_fusion_power_processing_and_conversion(self): + + # returns the tally with normalisation per pulse and conversion to joules + result = opp.process_damage_energy_tally( + source_strength=1.3e6, tally=self.my_tally, required_units="joule / second" + ) + + assert len(result) == 2 + assert result[0].units == "joule / second" + assert result[1].units == "joule / second" + + def test_energy_pulse_processing_and_conversion(self): + + # returns the tally with normalisation per pulse and conversion to joules + result = opp.process_damage_energy_tally( + source_strength=1.3e6, + tally=self.my_tally, + required_units="joules / pulse", # joules or joule can be requested + ) + + assert len(result) == 2 + assert result[0].units == "joule / pulse" + assert result[1].units == "joule / pulse" + + def test_energy_per_atom(self): + """makes use of material,""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="MeV / simulated_particle / atom", + volume=5, + material=my_mat, + ) + + assert len(result) == 2 + assert result[0].units == "megaelectron_volt / atom / simulated_particle" + assert result[1].units == "megaelectron_volt / atom / simulated_particle" + + def test_energy_per_second_per_atom(self): + """makes use of material and volume to convert to per atom units""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="MeV / atom / second", + source_strength=1, + volume=5, + material=my_mat, + ) + + assert len(result) == 2 + assert result[0].units == "megaelectron_volt / atom / second" + assert result[1].units == "megaelectron_volt / atom / second" + + def test_displacements_per_simulated_particle(self): + """makes use of energy_per_displacement to convert to displacements units""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / simulated_particle", + energy_per_displacement=80, + ) + + assert len(result) == 2 + assert result[0].units == "displacements / simulated_particle" + assert result[1].units == "displacements / simulated_particle" + + def test_displacements_per_second(self): + """makes use of energy_per_displacement to get displacements per second""" + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / second", + energy_per_displacement=80, + source_strength=100, + ) + + assert len(result) == 2 + assert result[0].units == "displacements / second" + assert result[1].units == "displacements / second" + + def test_displacements_per_atom(self): + """makes use of material and volume to find damage energy per atom""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / atom", + energy_per_displacement=80, + volume=5, + material=my_mat, + ) + + assert len(result) == 2 + assert result[0].units == "displacements / atom" + assert result[1].units == "displacements / atom" + + def test_displacements_per_atom_per_second(self): + """makes use of material and volume to find damage energy per atom per second""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / atom / second", + energy_per_displacement=80, + source_strength=100, + volume=5, + material=my_mat, + ) + + assert len(result) == 2 + assert result[0].units == "displacements / atom / second" + assert result[1].units == "displacements / atom / second" + + def test_displacements_per_atom_per_second_vs_per_year(self): + """makes use of material and volume to find damage energy per atom and + makes comparision between a year and second dpa""" + + my_mat = openmc.Material() + my_mat.add_element("Fe", 1) + my_mat.set_density("g/cm3", 1) + + result_second = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / atom / second", + energy_per_displacement=80, + source_strength=100, + volume=5, + material=my_mat, + ) + + assert len(result_second) == 2 + assert result_second[0].units == "displacements / atom / second" + assert result_second[1].units == "displacements / atom / second" + + result_year = opp.process_damage_energy_tally( + tally=self.my_tally, + required_units="displacements / atom / year", + energy_per_displacement=80, + source_strength=100, + volume=5, + material=my_mat, + ) + + assert len(result_year) == 2 + assert result_year[0].units == "displacements / atom / year" + assert result_year[1].units == "displacements / atom / year" + + assert ( + result_year[0].magnitude.sum() + == 365.25 * 24 * 60 * 60 * result_second[0].magnitude.sum() + ) + assert ( + result_year[1].magnitude.sum() + == 365.25 * 24 * 60 * 60 * result_second[1].magnitude.sum() + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_cell_heating.py b/tests/test_cell_heating.py index ea7a709..0ff24a5 100644 --- a/tests/test_cell_heating.py +++ b/tests/test_cell_heating.py @@ -1,4 +1,3 @@ - import unittest import openmc_post_processor as opp @@ -16,8 +15,7 @@ def setUp(self): def test_cell_tally_heating_no_processing(self): # returns the tally with base units result = opp.process_tally( - tally=self.my_tally, - required_units='eV / simulated_particle' + tally=self.my_tally, required_units="eV / simulated_particle" ) assert len(result) == 2 diff --git a/tests/test_cell_spectra.py b/tests/test_cell_spectra.py index 4fc44c4..8ad28c0 100644 --- a/tests/test_cell_spectra.py +++ b/tests/test_cell_spectra.py @@ -4,6 +4,7 @@ import pytest import openmc + class TestUsage(unittest.TestCase): def setUp(self): diff --git a/tests/test_python_api.py b/tests/test_python_api.py index 54a6bc6..c724d7f 100644 --- a/tests/test_python_api.py +++ b/tests/test_python_api.py @@ -4,6 +4,7 @@ import pytest import openmc + class TestUsage(unittest.TestCase): def setUp(self):