Merge pull request #17 from fusion-energy/develop

fixed a bug with volume normalisation

examples/processing_cell_flux_tally.py

@@ -10,7 +10,10 @@
 # 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)
+result = opp.process_tally(
+    tally=my_tally,
+    required_units="centimeter / simulated_particle"
+)
 print(f"flux base units = {result[0].units}", end="\n\n")
 
 

examples/processing_cell_heating_tally.py

@@ -9,7 +9,10 @@
 
 
 # returns the tally with base units
-result = opp.process_tally(tally=my_tally)
+result = opp.process_tally(
+    tally=my_tally,
+    required_units='eV / simulated_particle'
+)
 print(f"heating base units = {result}", end="\n\n")
 
 

examples/processing_mutliple_cell_spectra_tally.py

@@ -1,34 +1,54 @@
 import openmc
 import openmc_post_processor as opp
-from spectrum_plotter import plot_spectrum  # a convenient plotting package
+from spectrum_plotter import plot_spectrum_from_tally  # a convenient plotting package
 
 # loads in the statepoint file containing tallies
 statepoint = openmc.StatePoint(filepath="statepoint.2.h5")
 
+# constructs a dictionary of tallies
 results = {}
-for tally_name in ["2_neutron_spectra", "3_neutron_spectra"]:
+results["2_neutron_spectra"] = statepoint.get_tally(name="2_neutron_spectra")
+results["3_neutron_spectra"] = statepoint.get_tally(name="3_neutron_spectra")
 
-    # gets one tally from the available tallies
-    my_tally = statepoint.get_tally(name=tally_name)
-
-    # returns the tally with normalisation per pulse
-    result = opp.process_spectra_tally(
-        tally=my_tally,
-        required_units="centimeter / pulse",
-        required_energy_units="MeV",
-        source_strength=1.3e6,
-    )
-    results[tally_name] = result
+# plots a graph of the results with the units as recorded in the tally
+plot = plot_spectrum_from_tally(
+    spectrum=results,
+    x_label="Energy [MeV]",
+    y_label="neutron flux [centimeters / simulated_particle]",
+    x_scale="log",
+    y_scale="log",
+    filename="combine_spectra_plot_1.html",
+    required_energy_units="eV",
+    plotting_package='plotly',
+    required_units="centimeters / simulated_particle",
+)
 
 
-# plots a graph of the results
-plot = plot_spectrum(
+# plots a graph of the results with the units normalized for source strength
+plot = plot_spectrum_from_tally(
     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.html",
-    plotting_package='plotly'
+    filename="combine_spectra_plot_2.html",
+    required_energy_units="MeV",
+    plotting_package='plotly',
+    required_units="centimeter / second",
+    source_strength=1.3e6,
+)
+
+# plots a graph of the results with the units normalized for source strength and volume
+plot = plot_spectrum_from_tally(
+    spectrum=results,
+    x_label="Energy [MeV]",
+    y_label="neutron flux [neutrons s^-1 cm^-2]",
+    x_scale="log",
+    y_scale="log",
+    filename="combine_spectra_plot_3.html",
+    required_energy_units="MeV",
+    plotting_package='plotly',
+    required_units="neutrons / second * cm ** -2",
+    source_strength=1.3e7,
+    volume=100
 )

openmc_post_processor/utils.py

@@ -11,8 +11,8 @@
 
 def process_spectra_tally(
     tally,
-    required_units: str = None,
-    required_energy_units: str = None,
+    required_units: str = "centimeters / simulated_particle",
+    required_energy_units: str = 'eV',
     source_strength: float = None,
     volume: float = None,
 ) -> tuple:
@@ -22,7 +22,8 @@ def process_spectra_tally(
     Args:
         tally: The openmc.Tally object which should be a spectra tally. With a
             score of flux or current and an EnergyFilter
-        required_units: A tuple of the units to convert the energy and tally into
+        required_units: A tuple of 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
@@ -53,35 +54,35 @@ def process_spectra_tally(
     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,
-            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)
-        energy_in_required_units = energy_base.to(required_energy_units)
+    scaled_tally_result = scale_tally(
+        tally,
+        tally_base,
+        base_units,
+        ureg[required_units],
+        source_strength,
+        volume,
+    )
+    tally_in_required_units = scaled_tally_result.to(required_units)
+
+    scaled_tally_std_dev = scale_tally(
+        tally,
+        tally_std_dev_base,
+        base_units,
+        ureg[required_units],
+        source_strength,
+        volume,
+    )
 
-        return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units
-    else:
-        return energy_base, tally_base, tally_std_dev_base
+    tally_std_dev_in_required_units = scaled_tally_std_dev.to(required_units)
+
+    energy_in_required_units = energy_base.to(required_energy_units)
+
+    return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units
 
 
 def process_dose_tally(
     tally,
-    required_units: str = None,
+    required_units: str = 'picosievert cm **2 / simulated_particle',
     source_strength: float = None,
     volume: float = None,
 ):
@@ -119,33 +120,32 @@ def process_dose_tally(
     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,
-            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
+    scaled_tally_result = scale_tally(
+        tally,
+        tally_result,
+        base_units,
+        ureg[required_units],
+        source_strength,
+        volume,
+    )
+    tally_in_required_units = scaled_tally_result.to(required_units)
+
+    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 tally_result, tally_std_dev_base
+    return tally_in_required_units, tally_std_dev_in_required_units
 
 
 def process_tally(
     tally,
-    required_units: str = None,
+    required_units: str,
     source_strength: float = None,
     volume: float = None,
 ):
@@ -190,32 +190,36 @@ def process_tally(
     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,
-            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
+    scaled_tally_result = scale_tally(
+        tally,
+        tally_result,
+        base_units,
+        ureg[required_units],
+        source_strength,
+        volume,
+    )
+    tally_in_required_units = scaled_tally_result.to(required_units)
+
+    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 tally_result, tally_std_dev_base
+    return tally_in_required_units, tally_std_dev_in_required_units
 
 
 def scale_tally(
-    tally, tally_result, base_units, required_units, source_strength, volume
+    tally,
+    tally_result,
+    base_units,
+    required_units,
+    source_strength: float,
+    volume: float
 ):
     time_diff = check_for_dimentionality_difference(
         base_units, required_units, "[time]"
@@ -255,14 +259,16 @@ def scale_tally(
     if length_diff != 0:
         print("length scaling needed")
         if volume:
-            volume = volume * ureg["1 / centimeter ** 3"]
+            volume = volume * ureg["centimeter ** 3"]
         else:
             # volume required but not provided so it is found from the mesh
-            volume = compute_volume_of_voxels(tally) * ureg["1 / centimeter ** 3"]
+            volume = compute_volume_of_voxels(tally) * ureg["centimeter ** 3"]
 
-        if length_diff == -3:
-            tally_result = tally_result / volume
         if length_diff == 3:
+            print("dividing by volume")
+            tally_result = tally_result / volume
+        if length_diff == -3:
+            print("multiplying by volume")
             tally_result = tally_result * volume
 
     return tally_result