From 2d13fc6785ceafacf87821a55225725d2761370d Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 14:40:16 +0000 Subject: [PATCH 01/16] added units for dpa --- openmc_post_processor/utils.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 735dac3..b3f1623 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -424,11 +424,15 @@ 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 From 568dd4428302ab2132f9fc948b5241ef17df8fc2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 14:49:50 +0000 Subject: [PATCH 02/16] added damage-energy example --- .../processing_cell_damage-energy_tally.py | 21 +++++++++++++++++++ tests/create_statepoint_file_for_testing.py | 11 ++++++++-- ...e_statepoint_file_for_testing_one_batch.py | 10 +++++++-- ...statepoint_file_for_testing_using_dagmc.py | 11 ++++++++-- 4 files changed, 47 insertions(+), 6 deletions(-) create mode 100644 examples/processing_cell_damage-energy_tally.py diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py new file mode 100644 index 0000000..3cd2290 --- /dev/null +++ b/examples/processing_cell_damage-energy_tally.py @@ -0,0 +1,21 @@ +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_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_tally(tally=my_tally, required_units="MeV / simulated_particle") +print(f"damage-energy scaled base units = {result}", end="\n\n") 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 From 03630dc1e942515c811b219f757d6014344c13e1 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 1 Nov 2021 15:01:23 +0000 Subject: [PATCH 03/16] [skip ci] Apply formatting changes --- ...processing_2d_mesh_effective_dose_tally.py | 3 +- .../processing_cell_damage-energy_tally.py | 9 +++--- examples/processing_cell_flux_tally.py | 3 +- examples/processing_cell_heating_tally.py | 5 +-- examples/processing_cell_spectra_tally.py | 4 +-- .../processing_mutliple_cell_spectra_tally.py | 8 ++--- openmc_post_processor/__init__.py | 2 +- openmc_post_processor/utils.py | 31 ++++++++++++------- setup.py | 25 +++++++-------- tests/test_cell_heating.py | 4 +-- tests/test_cell_spectra.py | 1 + tests/test_python_api.py | 1 + 12 files changed, 48 insertions(+), 48 deletions(-) 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 index 3cd2290..16913f5 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -9,13 +9,12 @@ # 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"damage-energy base units = {result}", end="\n\n") # returns the tally with scalled based units (MeV instead of eV) -result, std_dev = opp.process_tally(tally=my_tally, required_units="MeV / simulated_particle") +result, std_dev = opp.process_tally( + tally=my_tally, required_units="MeV / simulated_particle" +) print(f"damage-energy scaled base units = {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..50e2fd0 100644 --- a/openmc_post_processor/__init__.py +++ b/openmc_post_processor/__init__.py @@ -7,5 +7,5 @@ process_tally, process_dose_tally, process_spectra_tally, - scale_tally + scale_tally, ) diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index b3f1623..b32436f 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -12,7 +12,7 @@ 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 +38,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() @@ -76,8 +76,12 @@ def process_spectra_tally( 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 +90,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,7 +115,7 @@ 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) @@ -230,7 +234,7 @@ def scale_tally( base_units, required_units, source_strength: float, - volume: float + volume: float, ): time_diff = check_for_dimentionality_difference( base_units, required_units, "[time]" @@ -377,13 +381,14 @@ def get_cell_ids_from_tally_filters(tally): # # 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 @@ -429,9 +434,11 @@ def get_tally_units(tally): units = ureg.electron_volt / ureg.simulated_particle else: - 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") + 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 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/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): From 9c5f6ac0e21b45aa8a69d377c5eaf60908e199d3 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 18:26:10 +0000 Subject: [PATCH 04/16] added method of convertingdamage-energy to dpa --- .../processing_cell_damage-energy_tally.py | 53 ++++- openmc_post_processor/__init__.py | 1 + openmc_post_processor/neutronics_units.txt | 2 + openmc_post_processor/utils.py | 209 ++++++++++++++---- 4 files changed, 209 insertions(+), 56 deletions(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 3cd2290..714f5c0 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -8,14 +8,51 @@ my_tally = statepoint.get_tally(name="2_damage-energy") -# returns the tally with base units -result = opp.process_tally( - tally=my_tally, - required_units='eV / simulated_particle' -) -print(f"damage-energy base units = {result}", end="\n\n") +# # 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") # returns the tally with scalled based units (MeV instead of eV) -result, std_dev = opp.process_tally(tally=my_tally, required_units="MeV / simulated_particle") -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 +) +print(f"damage-energy = {result}", end="\n\n") diff --git a/openmc_post_processor/__init__.py b/openmc_post_processor/__init__.py index 1887241..1f60e48 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, + process_damage_energy_tally, scale_tally ) diff --git a/openmc_post_processor/neutronics_units.txt b/openmc_post_processor/neutronics_units.txt index 4e771af..175fb57 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] +displacement = [displacement] pulse = [pulse] diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index b3f1623..330b3b7 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -9,6 +9,93 @@ 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 = None, + 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"]) * 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, + ) + if recombination_fraction: + scaled_tally_result = scaled_tally_result * recombination_fraction + 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", @@ -58,7 +145,6 @@ def process_spectra_tally( scaled_tally_result = scale_tally( tally, tally_base, - base_units, ureg[required_units], source_strength, volume, @@ -70,7 +156,6 @@ def process_spectra_tally( scaled_tally_std_dev = scale_tally( tally, tally_std_dev_base, - base_units, ureg[required_units], source_strength, volume, @@ -117,6 +202,11 @@ def process_dose_tally( 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 +216,6 @@ def process_dose_tally( scaled_tally_result = scale_tally( tally, tally_result, - base_units, ureg[required_units], source_strength, volume, @@ -138,7 +227,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 +286,6 @@ def process_tally( scaled_tally_result = scale_tally( tally, tally_result, - base_units, ureg[required_units], source_strength, volume, @@ -210,7 +297,6 @@ def process_tally( scaled_tally_std_dev = scale_tally( tally, tally_std_dev_base, - base_units, ureg[required_units], source_strength, volume, @@ -224,16 +310,47 @@ def process_tally( return tally_in_required_units + 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, "[displacement]" + ) + + print(time_diff, mass_diff, length_diff, displacement_diff) + input() + + 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['displacement'] + 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)") @@ -248,8 +365,9 @@ def scale_tally( f"source_strength is required but currently set to {source_strength}" ) + 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 +375,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,23 +383,43 @@ 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_with_units = compute_volume_of_voxels(tally) * ureg["centimeter ** 3"] 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 @@ -361,29 +499,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 +508,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 @@ -438,6 +549,8 @@ def get_tally_units(tally): 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 + # print('comparing ', units_1) + # print('comparing ', units_2) + 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 From eba135b5de3b17e3e1e92cc4066641b682e84048 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 1 Nov 2021 18:26:34 +0000 Subject: [PATCH 05/16] [skip ci] Apply formatting changes --- ...processing_2d_mesh_effective_dose_tally.py | 3 +- .../processing_cell_damage-energy_tally.py | 6 +- examples/processing_cell_flux_tally.py | 3 +- examples/processing_cell_heating_tally.py | 5 +- examples/processing_cell_spectra_tally.py | 4 +- .../processing_mutliple_cell_spectra_tally.py | 8 +-- openmc_post_processor/__init__.py | 2 +- openmc_post_processor/utils.py | 63 +++++++++++-------- setup.py | 25 ++++---- tests/test_cell_heating.py | 4 +- tests/test_cell_spectra.py | 1 + tests/test_python_api.py | 1 + 12 files changed, 66 insertions(+), 59 deletions(-) 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 index 714f5c0..823b411 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -35,8 +35,8 @@ # returns the tally with scalled based units (MeV instead of eV) my_mat = openmc.Material() -my_mat.add_element('Fe', 1) -my_mat.set_density('g/cm3', 1) +my_mat.add_element("Fe", 1) +my_mat.set_density("g/cm3", 1) # result, std_dev = opp.process_damage_energy_tally( # tally=my_tally, @@ -53,6 +53,6 @@ source_strength=1, volume=5, material=my_mat, - energy_per_displacement=10 + energy_per_displacement=10, ) 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 1f60e48..002bf1b 100644 --- a/openmc_post_processor/__init__.py +++ b/openmc_post_processor/__init__.py @@ -8,5 +8,5 @@ process_dose_tally, process_spectra_tally, process_damage_energy_tally, - scale_tally + scale_tally, ) diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 330b3b7..35c1e1e 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -9,10 +9,9 @@ ureg.load_definitions(str(Path(__file__).parent / "neutronics_units.txt")) - def process_damage_energy_tally( tally, - required_units: str = 'eV / simulated_particle', + required_units: str = "eV / simulated_particle", source_strength: float = None, volume: float = None, energy_per_displacement: float = None, @@ -45,7 +44,7 @@ def process_damage_energy_tally( """ if check_for_energy_function_filter(tally): - raise ValueError('EnergyFunctionFilter found in a damage-energy tally') + raise ValueError("EnergyFunctionFilter found in a damage-energy tally") # checks for user provided base units base_units = get_tally_units(tally) @@ -57,12 +56,12 @@ def process_damage_energy_tally( tally_result = np.array(data_frame["mean"]) * base_units if material: - atomic_mass_in_g = material.average_molar_mass * 1.66054E-24 + 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) + print("number_of_atoms_per_cm3", number_of_atoms_per_cm3) scaled_tally_result = scale_tally( tally, @@ -99,7 +98,7 @@ def process_damage_energy_tally( 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: @@ -125,10 +124,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() @@ -161,8 +160,12 @@ def process_spectra_tally( 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: @@ -171,7 +174,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, ): @@ -196,7 +199,7 @@ 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) @@ -310,7 +313,6 @@ def process_tally( return tally_in_required_units - def scale_tally( tally, tally_result, @@ -318,7 +320,7 @@ def scale_tally( source_strength: float, volume: float, atoms: float = None, - energy_per_displacement: float = None + energy_per_displacement: float = None, ): # energy_per_displacement @@ -338,17 +340,23 @@ def scale_tally( print(time_diff, mass_diff, length_diff, displacement_diff) input() - if time_diff == -2 and mass_diff==1 and length_diff == 2 and displacement_diff == -1: + 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['displacement'] + energy_per_displacement = ( + energy_per_displacement * ureg.electron_volt / ureg["displacement"] + ) 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]" ) @@ -365,7 +373,6 @@ def scale_tally( f"source_strength is required but currently set to {source_strength}" ) - time_diff = check_for_dimentionality_difference( tally_result.units, required_units, "[pulse]" ) @@ -391,7 +398,9 @@ def scale_tally( volume_with_units = volume * ureg["centimeter ** 3"] else: # volume required but not provided so it is found from the mesh - volume_with_units = compute_volume_of_voxels(tally) * ureg["centimeter ** 3"] + volume_with_units = ( + compute_volume_of_voxels(tally) * ureg["centimeter ** 3"] + ) if length_diff == 3: print("dividing by volume") @@ -416,9 +425,11 @@ def scale_tally( 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") + 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 @@ -540,9 +551,11 @@ def get_tally_units(tally): units = ureg.electron_volt / ureg.simulated_particle else: - 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") + 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 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/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): From 345df43cdccdb06a5d00d03b9212a41652681042 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 18:40:21 +0000 Subject: [PATCH 06/16] added recombination fraction --- examples/processing_cell_damage-energy_tally.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 714f5c0..2c0d2df 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -53,6 +53,7 @@ source_strength=1, volume=5, material=my_mat, - energy_per_displacement=10 + energy_per_displacement=10, + recombination_fraction=0.1 ) print(f"damage-energy = {result}", end="\n\n") From 3c5f1c9d6ef8d1f3c06a1c0684fe94a1a1a37471 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 18:41:10 +0000 Subject: [PATCH 07/16] added recombination factor --- openmc_post_processor/utils.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 330b3b7..56dad3c 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -16,7 +16,7 @@ def process_damage_energy_tally( source_strength: float = None, volume: float = None, energy_per_displacement: float = None, - recombination_fraction: float = None, + recombination_fraction: float = 0, material: float = None, ): """Processes a damage-energy tally converting the tally with default units @@ -54,7 +54,18 @@ def process_damage_energy_tally( print(f"tally {tally.name} base units {base_units}") - tally_result = np.array(data_frame["mean"]) * 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. - recombination_fraction) + + + tally_result = tally_result * base_units if material: atomic_mass_in_g = material.average_molar_mass * 1.66054E-24 @@ -73,8 +84,7 @@ def process_damage_energy_tally( number_of_atoms_per_cm3, energy_per_displacement, ) - if recombination_fraction: - scaled_tally_result = scaled_tally_result * recombination_fraction + tally_in_required_units = scaled_tally_result.to(required_units) if "std. dev." in data_frame.columns.to_list(): @@ -88,8 +98,8 @@ def process_damage_energy_tally( number_of_atoms_per_cm3, energy_per_displacement, ) - if recombination_fraction: - scaled_tally_std_dev = scaled_tally_std_dev * recombination_fraction + # 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: @@ -335,9 +345,6 @@ def scale_tally( tally_result.units, required_units, "[displacement]" ) - print(time_diff, mass_diff, length_diff, displacement_diff) - input() - 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: From 8c5af86616ff6e56515121c1ba9ebfdec8d867d2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 1 Nov 2021 18:58:40 +0000 Subject: [PATCH 08/16] [skip ci] merge conflict --- .../processing_cell_damage-energy_tally.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 3b1dbda..0ca179a 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -47,16 +47,24 @@ # ) # 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", - source_strength=1, + required_units="displacements per pulse", + source_strength=0.5, volume=5, material=my_mat, energy_per_displacement=10, -<<<<<<< HEAD recombination_fraction=0.1 -======= ->>>>>>> eba135b5de3b17e3e1e92cc4066641b682e84048 ) print(f"damage-energy = {result}", end="\n\n") From 3932ba32676db56eade5856cb73ba015af36509b Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 2 Nov 2021 09:34:19 +0000 Subject: [PATCH 09/16] [skip ci] Apply formatting changes --- examples/processing_cell_damage-energy_tally.py | 2 +- openmc_post_processor/utils.py | 15 +++++++++------ 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 0ca179a..7481164 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -65,6 +65,6 @@ volume=5, material=my_mat, energy_per_displacement=10, - recombination_fraction=0.1 + recombination_fraction=0.1, ) print(f"damage-energy = {result}", end="\n\n") diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 962e949..b34f2c7 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -53,17 +53,20 @@ def process_damage_energy_tally( print(f"tally {tally.name} base units {base_units}") - tally_result = np.array(data_frame["mean"]) - + 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}") + 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}") + raise ValueError( + f"recombination_fraction can't be larger than 1. recombination_fraction is {recombination_fraction}" + ) - tally_result = tally_result * (1. - recombination_fraction) + tally_result = tally_result * (1.0 - recombination_fraction) - tally_result = tally_result * base_units if material: From a03830d9f7227f7d7cc6beeec6620c4cfef6bcc2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Nov 2021 16:20:23 +0000 Subject: [PATCH 10/16] added tests for dpa --- .../processing_cell_damage-energy_tally.py | 75 +++++++++---------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 7481164..1cba112 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -8,55 +8,54 @@ 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 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 / 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") +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") -# returns the tally with scalled based units (MeV instead of eV) 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="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", + 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, From 837f2c1360718215529235820592893e81989528 Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 2 Nov 2021 16:20:41 +0000 Subject: [PATCH 11/16] [skip ci] Apply formatting changes --- examples/processing_cell_damage-energy_tally.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/processing_cell_damage-energy_tally.py b/examples/processing_cell_damage-energy_tally.py index 1cba112..e0fede1 100644 --- a/examples/processing_cell_damage-energy_tally.py +++ b/examples/processing_cell_damage-energy_tally.py @@ -10,8 +10,7 @@ # returns the tally with base units result = opp.process_damage_energy_tally( - tally=my_tally, - required_units='eV / simulated_particle' + tally=my_tally, required_units="eV / simulated_particle" ) print(f"damage-energy base units = {result}", end="\n\n") @@ -53,7 +52,7 @@ volume=5, material=my_mat, energy_per_displacement=10, - recombination_fraction=0.1 + recombination_fraction=0.1, ) print(f"damage-energy = {result}", end="\n\n") From 87b5afaa5cf18e1e16412106a022ecaf3dceb73c Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Nov 2021 16:22:47 +0000 Subject: [PATCH 12/16] added tests for dpa --- openmc_post_processor/neutronics_units.txt | 2 +- openmc_post_processor/utils.py | 39 +++- tests/test_cell_damage_energy.py | 257 +++++++++++++++++++++ 3 files changed, 289 insertions(+), 9 deletions(-) create mode 100644 tests/test_cell_damage_energy.py diff --git a/openmc_post_processor/neutronics_units.txt b/openmc_post_processor/neutronics_units.txt index 175fb57..8388bde 100644 --- a/openmc_post_processor/neutronics_units.txt +++ b/openmc_post_processor/neutronics_units.txt @@ -6,5 +6,5 @@ neutron = [] photon = [] simulated_particle = [] atom = [atom] -displacement = [displacement] +displacements = [displacements] pulse = [pulse] diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index b34f2c7..5b74d3c 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -347,7 +347,7 @@ def scale_tally( tally_result.units, required_units, "[mass]" ) displacement_diff = check_for_dimentionality_difference( - tally_result.units, required_units, "[displacement]" + tally_result.units, required_units, "[displacements]" ) if ( @@ -360,7 +360,7 @@ def scale_tally( print("energy per displacement_diff scaling needed (eV)") if energy_per_displacement: energy_per_displacement = ( - energy_per_displacement * ureg.electron_volt / ureg["displacement"] + energy_per_displacement * ureg.electron_volt / ureg["displacements"] ) tally_result = tally_result / energy_per_displacement else: @@ -368,6 +368,20 @@ def scale_tally( 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]" ) @@ -409,9 +423,16 @@ def scale_tally( volume_with_units = volume * ureg["centimeter ** 3"] else: # volume required but not provided so it is found from the mesh - volume_with_units = ( - 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") @@ -446,8 +467,9 @@ def scale_tally( 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] @@ -455,7 +477,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: 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() From de637fd66b2c673e03e26cfc3be5d695f07cc5f6 Mon Sep 17 00:00:00 2001 From: shimwell Date: Tue, 2 Nov 2021 16:23:56 +0000 Subject: [PATCH 13/16] [skip ci] Apply formatting changes --- openmc_post_processor/utils.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/openmc_post_processor/utils.py b/openmc_post_processor/utils.py index 5b74d3c..49bfa6c 100644 --- a/openmc_post_processor/utils.py +++ b/openmc_post_processor/utils.py @@ -375,7 +375,9 @@ def scale_tally( 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"] + energy_per_displacement = ( + energy_per_displacement * ureg.electron_volt / ureg["displacements"] + ) tally_result = tally_result / energy_per_displacement else: raise ValueError( @@ -428,10 +430,12 @@ def scale_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') + 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: From 27196f8d5b2292daa48c6eb2bcb57042f9cf11d8 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Nov 2021 16:49:13 +0000 Subject: [PATCH 14/16] added extra example and h5py --- .github/workflows/ci_with_install.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index 3074001..624b812 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -27,9 +27,10 @@ jobs: - name: install package run: | + python setup.py install + pip install h5py 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 run: | @@ -45,6 +46,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 From 80417eac1db093d2352c6fffa9323ecd6fcade10 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Nov 2021 16:56:47 +0000 Subject: [PATCH 15/16] trying conda update to fix h5py mismatch --- .github/workflows/ci_with_install.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index 624b812..9345b95 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -27,10 +27,10 @@ jobs: - name: install package run: | - python setup.py install - pip install h5py + conda update 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 run: | From f9f78c5e5ca2f7147472ab71843088ec0bc736c2 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 2 Nov 2021 17:06:48 +0000 Subject: [PATCH 16/16] removed conda update --- .github/workflows/ci_with_install.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.github/workflows/ci_with_install.yml b/.github/workflows/ci_with_install.yml index 9345b95..374a9a5 100644 --- a/.github/workflows/ci_with_install.yml +++ b/.github/workflows/ci_with_install.yml @@ -27,9 +27,7 @@ jobs: - name: install package run: | - conda update 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