Merge pull request #19 from fusion-energy/develop

allowing statepoints with a single batch to be processed

.github/workflows/ci_with_install.yml

@@ -38,6 +38,7 @@ jobs:
       - name: Run test_utils
         run: |
           python tests/create_statepoint_file_for_testing.py
+          python tests/create_statepoint_file_for_testing_one_batch.py
           pytest tests/ -v --cov=openmc_post_processor --cov-append --cov-report term --cov-report xml
 
       - name: Run examples

examples/processing_mutliple_cell_spectra_tally.py

@@ -1,5 +1,4 @@
 import openmc
-import openmc_post_processor as opp
 from spectrum_plotter import plot_spectrum_from_tally  # a convenient plotting package
 
 # loads in the statepoint file containing tallies
@@ -10,10 +9,11 @@
 results["2_neutron_spectra"] = statepoint.get_tally(name="2_neutron_spectra")
 results["3_neutron_spectra"] = statepoint.get_tally(name="3_neutron_spectra")
 
+
 # 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]",
+    x_label="Energy [eV]",
     y_label="neutron flux [centimeters / simulated_particle]",
     x_scale="log",
     y_scale="log",

openmc_post_processor/utils.py

@@ -50,10 +50,11 @@ def process_spectra_tally(
     print(f"tally {tally.name} base units {base_units}")
 
     # numpy array is needed as a pandas series can't have units
-    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
+    energy_in_required_units = energy_base.to(required_energy_units)
 
+    tally_base = np.array(data_frame["mean"]) * base_units
     scaled_tally_result = scale_tally(
         tally,
         tally_base,
@@ -64,20 +65,23 @@ def process_spectra_tally(
     )
     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)
+    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,
+            base_units,
+            ureg[required_units],
+            source_strength,
+            volume,
+        )
+        tally_std_dev_in_required_units = scaled_tally_std_dev.to(required_units)
+
+        return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units
 
-    energy_in_required_units = energy_base.to(required_energy_units)
+    else:
 
-    return energy_in_required_units, tally_in_required_units, tally_std_dev_in_required_units
+        return energy_in_required_units, tally_in_required_units
 
 
 def process_dose_tally(
@@ -118,7 +122,6 @@ def process_dose_tally(
     print(f"tally {tally.name} base units {base_units}")
 
     tally_result = np.array(data_frame["mean"]) * base_units
-    tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units
 
     scaled_tally_result = scale_tally(
         tally,
@@ -130,17 +133,20 @@ def process_dose_tally(
     )
     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_in_required_units, tally_std_dev_in_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,
+            base_units,
+            ureg[required_units],
+            source_strength,
+            volume,
+        )
+        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_tally(
@@ -188,7 +194,6 @@ def process_tally(
     print(f"tally {tally.name} base units {base_units}")
 
     tally_result = np.array(data_frame["mean"]) * base_units
-    tally_std_dev_base = np.array(data_frame["std. dev."]) * base_units
 
     scaled_tally_result = scale_tally(
         tally,
@@ -200,17 +205,23 @@ def process_tally(
     )
     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)
+    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,
+            base_units,
+            ureg[required_units],
+            source_strength,
+            volume,
+        )
+        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, tally_std_dev_in_required_units
+        return tally_in_required_units
 
 
 def scale_tally(

tests/create_statepoint_file_for_testing_one_batch.py

@@ -0,0 +1,166 @@
+# This minimal example makes a 3D volume and exports the shape to a stp file
+# A surrounding volume called a graveyard is needed for neutronics simulations
+
+import openmc
+import openmc_dagmc_wrapper as odw
+import openmc_plasma_source as ops
+import openmc_data_downloader as odd
+
+
+# MATERIALS
+breeder_material = openmc.Material(1, "PbLi")  # Pb84.2Li15.8
+breeder_material.add_element("Pb", 84.2, percent_type="ao")
+breeder_material.add_element(
+    "Li",
+    15.8,
+    percent_type="ao",
+    enrichment=50.0,
+    enrichment_target="Li6",
+    enrichment_type="ao",
+)  # 50% enriched
+breeder_material.set_density("atom/b-cm", 3.2720171e-2)  # around 11 g/cm3
+
+iron = openmc.Material(name="iron")
+iron.set_density("g/cm3", 7.75)
+iron.add_element("Fe", 0.95, percent_type="wo")
+
+materials = openmc.Materials([breeder_material, iron])
+
+odd.just_in_time_library_generator(libraries="TENDL-2019", materials=materials)
+
+# GEOMETRY
+
+# surfaces
+vessel_inner = openmc.Sphere(r=500)
+first_wall_outer_surface = openmc.Sphere(r=510)
+breeder_blanket_outer_surface = openmc.Sphere(r=610, boundary_type="vacuum")
+
+
+# cells
+inner_vessel_region = -vessel_inner
+inner_vessel_cell = openmc.Cell(region=inner_vessel_region)
+
+first_wall_region = -first_wall_outer_surface & +vessel_inner
+first_wall_cell = openmc.Cell(region=first_wall_region)
+first_wall_cell.fill = iron
+
+breeder_blanket_region = +first_wall_outer_surface & -breeder_blanket_outer_surface
+breeder_blanket_cell = openmc.Cell(region=breeder_blanket_region)
+breeder_blanket_cell.fill = breeder_material
+
+universe = openmc.Universe(
+    cells=[inner_vessel_cell, first_wall_cell, breeder_blanket_cell]
+)
+geometry = openmc.Geometry(universe)
+
+
+tally1 = odw.CellTally(
+    tally_type="flux",
+    target=2,
+)
+
+tally2 = odw.CellTally(
+    tally_type="neutron_spectra",
+    target=2,
+)
+
+tally3 = odw.CellTally(
+    tally_type="TBR",
+    target=2,
+)
+
+tally4 = odw.CellTally(
+    tally_type="heating",
+    target=2,
+)
+
+tally5 = odw.CellTally(
+    tally_type="neutron_effective_dose",
+    target=2,
+)
+
+tally6 = odw.CellTally(
+    tally_type="photon_effective_dose",
+    target=2,
+)
+
+tally7 = odw.CellTally(
+    tally_type="neutron_fast_flux",
+    target=2,
+)
+
+tally8 = odw.CellTally(
+    tally_type="photon_fast_flux",
+    target=2,
+)
+
+tally9 = odw.CellTally(
+    tally_type="neutron_spectra",
+    target=2,
+)
+
+tally10 = odw.CellTally(
+    tally_type="neutron_spectra",
+    target=2,
+)
+
+tally11 = odw.CellTally(
+    tally_type="neutron_spectra",
+    target=3,
+)
+
+tally12 = odw.CellTally(
+    tally_type="photon_spectra",
+    target=2,
+)
+
+tally13 = odw.CellTally(
+    tally_type="(n,total)",
+    target=2,
+)
+
+tally14 = odw.MeshTally2D(
+    tally_type="neutron_effective_dose",
+    plane="xy",
+    mesh_resolution=(10, 5),
+    bounding_box=[(-100, -100, 0), (100, 100, 1)],
+)
+
+tally15 = odw.MeshTally3D(
+    mesh_resolution=(100, 100, 100),
+    bounding_box=[(-100, -100, 0), (100, 100, 1)],
+    tally_type="neutron_effective_dose",
+)
+
+tallies = openmc.Tallies(
+    [
+        tally1,
+        tally2,
+        tally3,
+        tally4,
+        tally5,
+        tally6,
+        tally7,
+        tally8,
+        tally9,
+        tally10,
+        tally11,
+        tally12,
+        tally13,
+        tally14,
+        tally15,
+    ]
+)
+
+settings = odw.FusionSettings()
+settings.batches = 1
+settings.particles = 1000000
+# assigns a ring source of DT energy neutrons to the source using the
+# openmc_plasma_source package
+settings.source = ops.FusionPointSource()
+
+
+my_model = openmc.model.Model(
+    materials=materials, geometry=geometry, settings=settings, tallies=tallies
+)
+statepoint_file = my_model.run()