Skip to content

Commit

Permalink
Merge pull request #16 from fusion-energy/develop
Browse files Browse the repository at this point in the history
no longer reshaping values
shimwell authored Oct 27, 2021

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
2 parents ad3dd9e + c282c06 commit 510155a
Showing 6 changed files with 149 additions and 122 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -136,3 +136,4 @@ dmypy.json
*.xml
*.vscode/
*.png
*.html
10 changes: 4 additions & 6 deletions examples/processing_cell_flux_tally.py
Original file line number Diff line number Diff line change
@@ -7,20 +7,18 @@
# gets one tally from the available tallies
my_tally = statepoint.get_tally(name="2_flux")

print("default openmc tally result in base units", my_tally, end="\n\n")

print(opp.find_source_strength(fusion_energy_per_second_or_per_pulse=1.3e6))
# print(opp.find_source_strength(fusion_energy_per_second_or_per_pulse=1.3e6))

# returns the tally with base units
result = opp.process_tally(tally=my_tally)
print(f"flux base units = {result}", end="\n\n")
print(f"flux base units = {result[0].units}", end="\n\n")


# returns the tally with scalled based units (m instead of cm)
result = opp.process_tally(
tally=my_tally, required_units="meter / simulated_particle"
tally=my_tally, required_units="centimeter * particle / simulated_particle"
)
print(f"flux scaled base units = {result}", end="\n\n")
print(f"flux scaled base units = {result[0].units}", end="\n\n")


# returns the tally with normalisation per pulse
19 changes: 16 additions & 3 deletions examples/processing_cell_spectra_tally.py
Original file line number Diff line number Diff line change
@@ -15,18 +15,31 @@
)
print(f"spectra with base units = {result}", end="\n\n")

# returns the tally with base units
result = opp.process_spectra_tally(
tally=my_tally,
required_units="centimeter",
required_energy_units="MeV",
)
print(f"spectra with base units = {result}", end="\n\n")


# returns the tally with normalisation per pulse
result = opp.process_spectra_tally(
tally=my_tally, required_units=["eV", "centimeter / pulse"], source_strength=1.3e6
tally=my_tally,
required_units="centimeter / pulse",
required_energy_units="eV",
source_strength=1.3e6
)
print(f"spectra per pulse = {result}", end="\n\n")


# returns the tally scalled and normalisation for source strength
print(f"spectra per second = {result}", end="\n\n")
result = opp.process_spectra_tally(
tally=my_tally, required_units=["MeV", "centimeter / second"], source_strength=1e9
tally=my_tally,
required_units="centimeter / second",
required_energy_units="MeV",
source_strength=1e9
)
print(f"spectra per pulse = {result}", end="\n\n")

13 changes: 8 additions & 5 deletions examples/processing_mutliple_cell_spectra_tally.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,9 @@
import openmc
import openmc_post_processor as opp
from spectrum_plotter import plot_spectrum # a convenient plotting package

# loads in the statepoint file containing tallies
statepoint = opp.StatePoint(filepath="statepoint.2.h5")
statepoint = openmc.StatePoint(filepath="statepoint.2.h5")

results = {}
for tally_name in ["2_neutron_spectra", "3_neutron_spectra"]:
@@ -11,21 +12,23 @@
my_tally = statepoint.get_tally(name=tally_name)

# returns the tally with normalisation per pulse
result = statepoint.process_tally(
result = opp.process_spectra_tally(
tally=my_tally,
required_units=["MeV", "centimeter / pulse"],
required_units="centimeter / pulse",
required_energy_units="MeV",
source_strength=1.3e6,
)
results[tally_name] = result


# plots a graph of the results
plot_spectrum(
plot = plot_spectrum(
spectrum=results,
x_label="Energy [MeV]",
y_label="neutron flux [centimeter / second]",
x_scale="log",
y_scale="log",
# trim_zeros=False,
filename="combine_spectra_plot.png",
filename="combine_spectra_plot.html",
plotting_package='plotly'
)
2 changes: 0 additions & 2 deletions openmc_post_processor/__init__.py
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
from .utils import (
find_fusion_energy_per_reaction,
get_tally_units,
get_tally_units_dose,
get_tally_units_spectra,
check_for_dimentionality_difference,
find_source_strength,
compute_volume_of_voxels,
226 changes: 120 additions & 106 deletions openmc_post_processor/utils.py
Original file line number Diff line number Diff line change
@@ -11,7 +11,8 @@

def process_spectra_tally(
tally,
required_units: Tuple[str, str] = None,
required_units: str = None,
required_energy_units: str = None,
source_strength: float = None,
volume: float = None,
) -> tuple:
@@ -35,37 +36,47 @@ def process_spectra_tally(
Tuple of spectra energies and tally results
"""

check_for_unit_modifying_filters(tally)
if not check_for_energy_filter(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')

data_frame = tally.get_pandas_dataframe()

# checks for user provided base units
base_units = get_tally_units_spectra(tally)
base_units = get_tally_units(tally)
print(f"tally {tally.name} base units {base_units}")

for filter in tally.filters:
if isinstance(filter, openmc.filter.EnergyFilter):
energy_base = filter.values * base_units[0]
# skip other filters and contine or error if not found?

# numpy array is needed as a pandas series can't have units
tally_base = np.array(data_frame["mean"]) * base_units[1]
tally_base = np.array(data_frame["mean"]) * base_units
tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units
energy_base = np.array(data_frame["energy low [eV]"]) * ureg.electron_volt

if required_units:
scaled_tally_result = scale_tally(
tally,
tally_base,
base_units[1],
ureg[required_units[1]],
base_units,
ureg[required_units],
source_strength,
volume,
)
scaled_tally_std_dev = scale_tally(
tally,
tally_std_dev_base,
base_units,
ureg[required_units],
source_strength,
volume,
)
tally_in_required_units = scaled_tally_result.to(required_units[1])
energy_in_required_units = energy_base.to(required_units[0])
tally_std_dev_in_required_units = scaled_tally_std_dev.to(required_units)
tally_in_required_units = scaled_tally_result.to(required_units)
energy_in_required_units = energy_base.to(required_energy_units)

return energy_in_required_units, tally_in_required_units
return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units
else:
return energy_base, tally_base
return energy_base, tally_base, tally_std_dev_base


def process_dose_tally(
@@ -94,38 +105,42 @@ def process_dose_tally(
The dose tally result in the required units
"""

data_frame = tally.get_pandas_dataframe()
if not check_for_energy_function_filter(tally):
raise ValueError('EnergyFunctionFilter was not found in dose tally')

# checks for user provided base units
base_units = get_tally_units_dose(tally)
base_units = get_tally_units(tally)
base_units = ureg.picosievert * ureg.centimeter * base_units

data_frame = tally.get_pandas_dataframe()

print(f"tally {tally.name} base units {base_units}")

if len(data_frame["mean"]) == 1:
# just a single number in the tally result
tally_result = data_frame["mean"].sum() * base_units[0]
else:
# more than one number, a mesh tally or several cell tallies
tally_filter = tally.find_filter(filter_type=openmc.MeshFilter)
shape = tally_filter.mesh.dimension.tolist()
if 1 in shape:
# 2d mesh
shape.remove(1)
tally_result = np.array(data_frame["mean"]) * base_units[0]
tally_result = tally_result.reshape(shape)
tally_result = np.array(data_frame["mean"]) * base_units
tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units

if required_units:
scaled_tally_result = scale_tally(
tally,
tally_result,
base_units[0],
base_units,
ureg[required_units],
source_strength,
volume,
)
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)
tally_in_required_units = scaled_tally_result.to(required_units)
return tally_in_required_units
return tally_in_required_units, tally_std_dev_in_required_units

return tally_result
return tally_result, tally_std_dev_base


def process_tally(
@@ -152,42 +167,51 @@ def process_tally(
Returns:
The dose tally result in the required units
"""
data_frame = tally.get_pandas_dataframe()

check_for_unit_modifying_filters(tally)
if check_for_energy_function_filter(tally):
msg = (
"An EnergyFunctionFilter was found in the tally. This "
"modifies the tally units and the base units of the"
"EnergyFunctionFilter are not known to OpenMC. Therefore "
"the units of this tally can not be found. If you have "
"applied dose coefficients to an EnergyFunctionFilter "
"the units of these are known and yo can use the "
"get_tally_units_dose function instead of the "
"get_tally_units"
)
raise ValueError(msg)

data_frame = tally.get_pandas_dataframe()

# checks for user provided base units
base_units = get_tally_units(tally)

print(f"tally {tally.name} base units {base_units}")

# there might be more than one based unit entry if spectra has been tallied
if len(base_units) == 1:
if len(data_frame["mean"]) == 1:
# just a single number in the tally result
tally_result = data_frame["mean"].sum() * base_units[0]
else:
# more than one number, a mesh tally
tally_filter = tally.find_filter(filter_type=openmc.MeshFilter)
shape = tally_filter.mesh.dimension.tolist()
if 1 in shape:
# 2d mesh
shape.remove(1)
tally_result = np.array(data_frame["mean"]) * base_units[0]
tally_result = tally_result.reshape(shape)

if required_units:
scaled_tally_result = scale_tally(
tally,
tally_result,
base_units[0],
ureg[required_units],
source_strength,
volume,
)
tally_in_required_units = scaled_tally_result.to(required_units)
return tally_in_required_units
tally_result = np.array(data_frame["mean"]) * base_units
tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units

return tally_result
if required_units:
scaled_tally_result = scale_tally(
tally,
tally_result,
base_units,
ureg[required_units],
source_strength,
volume,
)
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)
tally_in_required_units = scaled_tally_result.to(required_units)
return tally_in_required_units, tally_std_dev_in_required_units

return tally_result, tally_std_dev_base


def scale_tally(
@@ -320,78 +344,68 @@ def get_cell_ids_from_tally_filters(tally):
return cell_ids


def get_tally_units_spectra(tally):
# def get_tally_units_dose(tally):

units = get_tally_units(tally)

# check it is a spectra tally by looking for a openmc.filter.EnergyFilter
for filter in tally.filters:
if isinstance(filter, openmc.filter.EnergyFilter):
# spectra tally has units for the energy as well as the flux
units = [ureg.electron_volt, units[0]]
return units
# # An EnergyFunctionFilter is exspeted in dose tallies
# units = get_tally_units(tally)

raise ValueError(
"units for spectra tally can't be found, an EnergyFilter was not present"
)
# 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

def get_tally_units_dose(tally):
# raise ValueError(
# "units for dose tally can't be found, an EnergyFunctionFilter was not present"
# )

# An EnergyFunctionFilter is exspeted in dose tallies
units = get_tally_units(tally)
def check_for_energy_filter(tally):

# check it is a dose tally by looking for a openmc.filter.EnergyFunctionFilter
# check it is a spectra tally by looking for a openmc.filter.EnergyFilter
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
units = [ureg.picosievert * ureg.centimeter * units[0]]
return units

raise ValueError(
"units for dose tally can't be found, an EnergyFunctionFilter was not present"
)
if isinstance(filter, openmc.filter.EnergyFilter):
# spectra tally has units for the energy as well as the flux
return True
return False

# raise ValueError(
# "units for spectra tally can't be found, an EnergyFilter was not present"
# )


def check_for_unit_modifying_filters(tally):
def check_for_energy_function_filter(tally):
# check for EnergyFunctionFilter which modify the units of the tally
for filter in tally.filters:
if isinstance(filter, openmc.filter.EnergyFunctionFilter):
msg = (
"An EnergyFunctionFilter was found in the tally. This "
"modifies the tally units and the base units of the"
"EnergyFunctionFilter are not known to OpenMC. Therefore "
"the units of this tally can not be found. If you have "
"applied dose coefficients to an EnergyFunctionFilter "
"the units of these are known and yo can use the "
"get_tally_units_dose function instead of the "
"get_tally_units"
)
raise ValueError(msg)
return True
return False


def get_tally_units(tally):
""" """

if tally.scores == ["current"]:
units = get_particles_from_tally_filters(tally, ureg)
units = [units / (ureg.simulated_particle * ureg.centimeter ** 2)]
units = units / (ureg.simulated_particle * ureg.centimeter ** 2)

if tally.scores == ["flux"]:
print("score is flux")
# tally has units of particle-cm2 per simulated_particle
# https://openmc.discourse.group/t/normalizing-tally-to-get-flux-value/99/4
units = get_particles_from_tally_filters(tally, ureg)
units = [units * ureg.centimeter / ureg.simulated_particle]
units = units * ureg.centimeter / ureg.simulated_particle

elif tally.scores == ["heating"]:
# heating units are eV / simulated_particle
units = [ureg.electron_volt / ureg.simulated_particle]
units = ureg.electron_volt / ureg.simulated_particle

else:
# return [1 / ureg.simulated_particle]

0 comments on commit 510155a

Please sign in to comment.