Merge pull request #5 from fusion-energy/adding-phi-r-plot

Adding phi r plot

.gitignore

@@ -132,3 +132,4 @@ src/*/_version.py
 *.vtk
 *.h5
 *.out
+*.png

README.md

@@ -5,6 +5,11 @@
 Plots R-Z or Phi-R slices of OpenMC cylindrical mesh tallies
 
 
+![openmc cylinder PhiR mesh tally plot](https://user-images.githubusercontent.com/8583900/227660135-8ac2eb69-829c-4788-ae02-f9e4b1465aa6.png)
+![openmc cylinder RZ mesh tally plot](https://user-images.githubusercontent.com/8583900/227660142-99bcb264-57f8-47d6-b9ff-d211645fe185.png)
+
+
+
 ## Installation
 
 openmc_cylindrical_mesh_plotter is available on PyPi

examples/plot_phir_slice_point_source.py

@@ -0,0 +1,74 @@
+# this example creates a simple CylindricalMesh tally and performs an openmc
+# simulation to populate the tally. Slices of the resulting tally is then
+# plotted using the openmc_cylindrical_plotter in the Phi R axis.
+
+import openmc
+import numpy as np
+from math import pi
+import matplotlib.pyplot as plt
+import openmc_cylindrical_mesh_plotter  # adds slice_of_data method to CylindricalMesh
+from matplotlib import ticker
+
+mesh = openmc.CylindricalMesh()
+mesh.phi_grid = np.linspace(0, 2 * pi, 50)  # note the mesh is full 2pi circle
+mesh.r_grid = np.linspace(0, 10, 20)
+mesh.z_grid = np.linspace(0, 5, 4)
+
+tally = openmc.Tally(name="my_tally")
+mesh_filter = openmc.MeshFilter(mesh)
+tally.filters.append(mesh_filter)
+tally.scores.append("flux")
+tallies = openmc.Tallies([tally])
+
+outer_surface = openmc.Sphere(r=100, boundary_type="vacuum")
+cell = openmc.Cell(region=-outer_surface)
+
+material = openmc.Material()
+material.add_nuclide("Fe56", 1)
+material.set_density("g/cm3", 0.1)
+my_materials = openmc.Materials([material])
+
+universe = openmc.Universe(cells=[cell])
+my_geometry = openmc.Geometry(universe)
+
+my_source = openmc.Source()
+
+# this makes a point source instead with
+my_source.space = openmc.stats.Point((0, -5, 0))
+# sets the direction to isotropic
+my_source.angle = openmc.stats.Isotropic()
+
+
+my_settings = openmc.Settings()
+# my_settings.inactive = 0
+my_settings.run_mode = "fixed source"
+my_settings.batches = 10
+my_settings.particles = 100000
+my_settings.source = my_source
+
+model = openmc.model.Model(my_geometry, my_materials, my_settings, tallies)
+sp_filename = model.run()
+
+statepoint = openmc.StatePoint(sp_filename)
+
+my_tally_result = statepoint.get_tally(name="my_tally")
+
+for slice_index in range(1, len(mesh.z_grid)):
+    theta, r, values = mesh.slice_of_data(
+        dataset=my_tally_result.mean.flatten(),
+        slice_index=slice_index,
+        axis="PhiR",
+        volume_normalization=False,
+    )
+
+    fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
+    im = ax.contourf(
+        theta, r, values, extent=(0, 100, 0, 50)
+    )  # , locator=ticker.LogLocator())
+
+    # sets the y axis limits to match the mesh limits
+    ax.set_ylim(mesh.r_grid[0], mesh.r_grid[-1])
+
+    plt.colorbar(im, label="Flux")
+
+    plt.show()

examples/plot_phir_slice_ring_source.py

@@ -0,0 +1,87 @@
+# this example creates a simple CylindricalMesh tally and performs an openmc
+# simulation to populate the tally. Slices of the resulting tally is then
+# plotted using the openmc_cylindrical_plotter
+
+import openmc
+import numpy as np
+from math import pi
+import matplotlib.pyplot as plt
+import openmc_cylindrical_mesh_plotter  # adds slice_of_data method to CylindricalMesh
+from matplotlib import ticker
+from math import pi
+
+mesh = openmc.CylindricalMesh()
+mesh.phi_grid = np.linspace(
+    0.0, 1.5 * pi, 10
+)  # note the mesh is 3/4 of a circle, not the full 2pi
+mesh.r_grid = np.linspace(0, 10, 20)
+mesh.z_grid = np.linspace(0, 5, 4)
+
+tally = openmc.Tally(name="my_tally")
+mesh_filter = openmc.MeshFilter(mesh)
+tally.filters.append(mesh_filter)
+tally.scores.append("flux")
+tallies = openmc.Tallies([tally])
+
+outer_surface = openmc.Sphere(r=100, boundary_type="vacuum")
+cell = openmc.Cell(region=-outer_surface)
+
+material = openmc.Material()
+material.add_nuclide("Fe56", 1)
+material.set_density("g/cm3", 0.1)
+my_materials = openmc.Materials([material])
+
+universe = openmc.Universe(cells=[cell])
+my_geometry = openmc.Geometry(universe)
+
+my_source = openmc.Source()
+
+# the distribution of radius is just a single value
+radius = openmc.stats.Discrete([5], [1])
+# the distribution of source z values is just a single value
+z_values = openmc.stats.Discrete([2.5], [1])
+# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi
+angle = openmc.stats.Uniform(a=0.0, b=pi)  # half the circle 0 to 180 degrees
+# this makes the ring source using the three distributions and a radius
+# could do a point source instead with
+# my_source.space = openmc.stats.Point((0,0,0))
+my_source.space = openmc.stats.CylindricalIndependent(
+    r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)
+)
+# sets the direction to isotropic
+my_source.angle = openmc.stats.Isotropic()
+
+
+my_settings = openmc.Settings()
+# my_settings.inactive = 0
+my_settings.run_mode = "fixed source"
+my_settings.batches = 10
+my_settings.particles = 100000
+my_settings.source = my_source
+
+model = openmc.model.Model(my_geometry, my_materials, my_settings, tallies)
+sp_filename = model.run()
+
+statepoint = openmc.StatePoint(sp_filename)
+
+my_tally_result = statepoint.get_tally(name="my_tally")
+
+for slice_index in range(1, len(mesh.z_grid)):
+    theta, r, values = mesh.slice_of_data(
+        dataset=my_tally_result.mean.flatten(),
+        slice_index=slice_index,
+        axis="PhiR",
+        volume_normalization=False,
+    )
+
+    fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
+    im = ax.contourf(
+        theta, r, values, extent=(0, 100, 0, 50)
+    )  # , locator=ticker.LogLocator())
+
+    # sets the y axis limits to match the mesh limits
+    ax.set_ylim(mesh.r_grid[0], mesh.r_grid[-1])
+
+    plt.colorbar(im, label="Flux")
+
+    plt.show()

examples/plot_rz_slices.py → examples/plot_rz_slices_point_source.py

@@ -61,9 +61,10 @@
 
 my_tally_result = statepoint.get_tally(name="my_tally")
 
-for slice_index in range(len(mesh.phi_grid) - 1):
+for slice_index in range(1, len(mesh.phi_grid)):
     data = mesh.slice_of_data(
-        dataset=my_tally_result.mean,
+        dataset=my_tally_result.mean.flatten(),
+        axis="RZ",
         # dataset=np.array(2*19*49*[1]), flat data for testing
         slice_index=slice_index,
         volume_normalization=False,
@@ -72,6 +73,8 @@
     x_label, y_label = mesh.get_axis_labels()
     plt.xlabel(x_label)
     plt.ylabel(y_label)
-    plt.imshow(data, extent=extent)
+    im = plt.imshow(data, extent=extent)
+
+    plt.colorbar(im, label="Flux")
 
     plt.show()

examples/plot_rz_slices_ring_source.py

@@ -0,0 +1,79 @@
+# this example creates a simple CylindricalMesh tally and performs an openmc
+# simulation to populate the tally. Slices of the resulting tally is then
+# plotted using the openmc_cylindrical_plotter in the R Z axis
+
+import openmc
+import numpy as np
+from math import pi
+import matplotlib.pyplot as plt
+import openmc_cylindrical_mesh_plotter  # adds slice_of_data method to CylindricalMesh
+
+mesh = openmc.CylindricalMesh()
+mesh.phi_grid = np.linspace(0.0, 2 * pi, 3)
+mesh.r_grid = np.linspace(0, 10, 50)
+mesh.z_grid = np.linspace(0, 8, 50)
+
+tally = openmc.Tally(name="my_tally")
+mesh_filter = openmc.MeshFilter(mesh)
+tally.filters.append(mesh_filter)
+tally.scores.append("flux")
+tallies = openmc.Tallies([tally])
+
+outer_surface = openmc.Sphere(r=100, boundary_type="vacuum")
+cell = openmc.Cell(region=-outer_surface)
+
+material = openmc.Material()
+material.add_nuclide("Fe56", 1)
+material.set_density("g/cm3", 0.1)
+my_materials = openmc.Materials([material])
+
+universe = openmc.Universe(cells=[cell])
+my_geometry = openmc.Geometry(universe)
+
+my_source = openmc.Source()
+
+# the distribution of radius is just a single value
+radius = openmc.stats.Discrete([5], [1])
+# the distribution of source z values is just a single value
+z_values = openmc.stats.Discrete([4], [1])
+# the distribution of source azimuthal angles values is a uniform distribution between 0 and 2 Pi
+angle = openmc.stats.Uniform(a=0.0, b=2 * 3.14159265359)
+# this makes the ring source using the three distributions and a radius
+# could do a point source instead with my_source.space = openmc.stats.Point((5,5., 5))
+my_source.space = openmc.stats.CylindricalIndependent(
+    r=radius, phi=angle, z=z_values, origin=(0.0, 0.0, 0.0)
+)
+# sets the direction to isotropic
+my_source.angle = openmc.stats.Isotropic()
+
+
+my_settings = openmc.Settings()
+# my_settings.inactive = 0
+my_settings.run_mode = "fixed source"
+my_settings.batches = 10
+my_settings.particles = 100000
+my_settings.source = my_source
+
+model = openmc.model.Model(my_geometry, my_materials, my_settings, tallies)
+sp_filename = model.run()
+
+statepoint = openmc.StatePoint(sp_filename)
+
+my_tally_result = statepoint.get_tally(name="my_tally")
+
+for slice_index in range(1, len(mesh.phi_grid)):
+    data = mesh.slice_of_data(
+        dataset=my_tally_result.mean.flatten(),
+        axis="RZ",
+        # dataset=np.array(2*19*49*[1]), flat data for testing
+        slice_index=slice_index,
+        volume_normalization=False,
+    )
+    extent = mesh.get_mpl_plot_extent()
+    x_label, y_label = mesh.get_axis_labels()
+    plt.title("Ring source with radius=5, z=4")
+    plt.xlabel(x_label)
+    plt.ylabel(y_label)
+    im = plt.imshow(data, extent=extent)
+    plt.colorbar(im, label="Flux")
+    plt.show()