From 0c40d99291015890137b9a1f20b03a2008663149 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 22 Mar 2023 17:10:01 +0000 Subject: [PATCH 1/9] example phi-r slice looks correct at first pass --- examples/plot_phir_slices.py | 90 ++++++++++++++++++++++++++++++------ 1 file changed, 77 insertions(+), 13 deletions(-) diff --git a/examples/plot_phir_slices.py b/examples/plot_phir_slices.py index 83b85a2..b68511b 100644 --- a/examples/plot_phir_slices.py +++ b/examples/plot_phir_slices.py @@ -1,31 +1,95 @@ -import numpy as np -import matplotlib.pyplot as plt -from math import pi +# 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 mesh = openmc.CylindricalMesh() mesh.phi_grid = np.linspace(0.0, 2 * pi, 10) -mesh.r_grid = np.linspace(0, 10, 4) -mesh.z_grid = np.linspace(0, 5, 5) +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=1 * 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((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() -sp_filename = "examples/statepoint.10.h5" statepoint = openmc.StatePoint(sp_filename) my_tally_result = statepoint.get_tally(name="my_tally") -actual = mesh.phi_grid +actual = np.linspace(mesh.phi_grid[0],mesh.phi_grid[-1], mesh.dimension[1]) +expected = np.linspace(mesh.r_grid[0],mesh.r_grid[-1], mesh.dimension[0]) # Using linspace so that the endpoint of 360 is included # actual = np.radians(np.linspace(0, 360, 20)) # expected = np.arange(0, 70, 10) -expected = mesh.r_grid - -# r, theta = np.meshgrid(expected, actual) +r, theta = np.meshgrid(expected, actual) +values=np.array([(len(mesh.r_grid)-1)*[12]]*(len(mesh.phi_grid)-1)) # # values = np.random.random((actual.size, expected.size)) -# fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) -# ax.contourf(theta, r, values) +for slice_index in range(len(mesh.z_grid)-1): + lower_index = int(slice_index*(len(mesh.phi_grid)-1)) + upper_index = int((slice_index+1)*(len(mesh.phi_grid)-1)) + + dataset=my_tally_result.mean + + + # values=dataset.flatten().reshape(-1,len(mesh.r_grid)-1,order='C')[:len(mesh.phi_grid)-1] + values=dataset.flatten().reshape(-1,len(mesh.r_grid)-1,order='A')[lower_index:upper_index] + # rot_val = np.rot90(values) + + fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) + from matplotlib import ticker + im = ax.contourf(theta, r, values)#, locator=ticker.LogLocator()) + + # ax.contour(theta, r, theta) + # ax.contourf(theta, r, data_slice) -# plt.show() + plt.colorbar(im) + plt.show() From 376f332d0f83f18695668c3efd48f3c4732ad430 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 22 Mar 2023 17:11:44 +0000 Subject: [PATCH 2/9] [skip ci] Apply formatting changes --- examples/plot_phir_slices.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/examples/plot_phir_slices.py b/examples/plot_phir_slices.py index b68511b..3e5d247 100644 --- a/examples/plot_phir_slices.py +++ b/examples/plot_phir_slices.py @@ -62,31 +62,33 @@ my_tally_result = statepoint.get_tally(name="my_tally") -actual = np.linspace(mesh.phi_grid[0],mesh.phi_grid[-1], mesh.dimension[1]) -expected = np.linspace(mesh.r_grid[0],mesh.r_grid[-1], mesh.dimension[0]) +actual = np.linspace(mesh.phi_grid[0], mesh.phi_grid[-1], mesh.dimension[1]) +expected = np.linspace(mesh.r_grid[0], mesh.r_grid[-1], mesh.dimension[0]) # Using linspace so that the endpoint of 360 is included # actual = np.radians(np.linspace(0, 360, 20)) # expected = np.arange(0, 70, 10) r, theta = np.meshgrid(expected, actual) -values=np.array([(len(mesh.r_grid)-1)*[12]]*(len(mesh.phi_grid)-1)) +values = np.array([(len(mesh.r_grid) - 1) * [12]] * (len(mesh.phi_grid) - 1)) # # values = np.random.random((actual.size, expected.size)) -for slice_index in range(len(mesh.z_grid)-1): - lower_index = int(slice_index*(len(mesh.phi_grid)-1)) - upper_index = int((slice_index+1)*(len(mesh.phi_grid)-1)) - - dataset=my_tally_result.mean +for slice_index in range(len(mesh.z_grid) - 1): + lower_index = int(slice_index * (len(mesh.phi_grid) - 1)) + upper_index = int((slice_index + 1) * (len(mesh.phi_grid) - 1)) + dataset = my_tally_result.mean # values=dataset.flatten().reshape(-1,len(mesh.r_grid)-1,order='C')[:len(mesh.phi_grid)-1] - values=dataset.flatten().reshape(-1,len(mesh.r_grid)-1,order='A')[lower_index:upper_index] + values = dataset.flatten().reshape(-1, len(mesh.r_grid) - 1, order="A")[ + lower_index:upper_index + ] # rot_val = np.rot90(values) - fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) + fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) from matplotlib import ticker - im = ax.contourf(theta, r, values)#, locator=ticker.LogLocator()) + + im = ax.contourf(theta, r, values) # , locator=ticker.LogLocator()) # ax.contour(theta, r, theta) # ax.contourf(theta, r, data_slice) From 49f86db22b78d27fc687518a8c77c2a660af4e09 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 23 Mar 2023 12:31:43 +0000 Subject: [PATCH 3/9] refactorign tests --- examples/plot_phir_slices.py | 36 ++--- src/openmc_cylindrical_mesh_plotter/core.py | 49 ++++++- tests/test_slice_of_data.py | 149 +++++++++++++------- 3 files changed, 154 insertions(+), 80 deletions(-) diff --git a/examples/plot_phir_slices.py b/examples/plot_phir_slices.py index 3e5d247..fdb941e 100644 --- a/examples/plot_phir_slices.py +++ b/examples/plot_phir_slices.py @@ -7,6 +7,7 @@ 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.0, 2 * pi, 10) @@ -62,36 +63,19 @@ my_tally_result = statepoint.get_tally(name="my_tally") -actual = np.linspace(mesh.phi_grid[0], mesh.phi_grid[-1], mesh.dimension[1]) -expected = np.linspace(mesh.r_grid[0], mesh.r_grid[-1], mesh.dimension[0]) - -# Using linspace so that the endpoint of 360 is included -# actual = np.radians(np.linspace(0, 360, 20)) -# expected = np.arange(0, 70, 10) - -r, theta = np.meshgrid(expected, actual) -values = np.array([(len(mesh.r_grid) - 1) * [12]] * (len(mesh.phi_grid) - 1)) -# # values = np.random.random((actual.size, expected.size)) - for slice_index in range(len(mesh.z_grid) - 1): - lower_index = int(slice_index * (len(mesh.phi_grid) - 1)) - upper_index = int((slice_index + 1) * (len(mesh.phi_grid) - 1)) - - dataset = my_tally_result.mean - - # values=dataset.flatten().reshape(-1,len(mesh.r_grid)-1,order='C')[:len(mesh.phi_grid)-1] - values = dataset.flatten().reshape(-1, len(mesh.r_grid) - 1, order="A")[ - lower_index:upper_index - ] - # rot_val = np.rot90(values) + theta, r, values = mesh.slice_of_data( + dataset=my_tally_result.mean, + slice_index=slice_index, + axis='Phi-R', + volume_normalization=False, + ) + # plt.xlabel(x_label) + # plt.ylabel(y_label) fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) - from matplotlib import ticker - im = ax.contourf(theta, r, values) # , locator=ticker.LogLocator()) - # ax.contour(theta, r, theta) - # ax.contourf(theta, r, data_slice) - plt.colorbar(im) + plt.show() diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index b489588..45fdabb 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -6,10 +6,57 @@ def slice_of_data( self, dataset: np.ndarray, - # axis='RZ' todo add more axis including top down view + axis: str, slice_index=0, volume_normalization: bool = True, ): + if axis == 'R-Z': + return slice_of_rz_data( + self, + dataset=dataset, + slice_index=slice_index, + volume_normalization=volume_normalization + ) + elif axis == 'Phi-R': + return slice_of_phir_data( + self, + dataset=dataset, + slice_index=slice_index, + volume_normalization=volume_normalization + ) + else: + raise ValueError(f'axis must be either "R-Z" or "Phi-R", not {axis}') + + +def slice_of_phir_data( + self, + dataset: np.ndarray, + slice_index=0, + volume_normalization: bool = True, +): + actual = np.linspace(self.phi_grid[0], self.phi_grid[-1], self.dimension[1]) + expected = np.linspace(self.r_grid[0], self.r_grid[-1], self.dimension[0]) + + r, theta = np.meshgrid(expected, actual) + + lower_index = int(slice_index * (len(self.phi_grid) - 1)) + upper_index = int((slice_index + 1) * (len(self.phi_grid) - 1)) + + # both order A and C appear to work + # values=dataset.flatten().reshape(-1,len(self.r_grid)-1,order='C')[:len(self.phi_grid)-1] + values = dataset.flatten().reshape(-1, len(self.r_grid) - 1, order="A")[ + lower_index:upper_index + ] + + return theta, r, values + +def slice_of_rz_data( + self, + dataset: np.ndarray, + slice_index=0, + volume_normalization: bool = True, +): + lower_index = int(slice_index * (len(self.r_grid) - 1)) upper_index = int((slice_index + 1) * (len(self.r_grid) - 1)) diff --git a/tests/test_slice_of_data.py b/tests/test_slice_of_data.py index 6eaf699..af032ee 100644 --- a/tests/test_slice_of_data.py +++ b/tests/test_slice_of_data.py @@ -5,64 +5,85 @@ import openmc_cylindrical_mesh_plotter # adds slice_of_data method to CylindricalMesh -def test_slice_of_data(): - mesh = openmc.CylindricalMesh() - mesh.phi_grid = np.linspace(0.0, 2 * pi, 10) - mesh.r_grid = np.linspace(0, 10, 4) - mesh.z_grid = np.linspace(0, 5, 5) - - 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=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") +mesh = openmc.CylindricalMesh() +mesh.phi_grid = np.linspace(0.0, 2 * pi, 10) +mesh.r_grid = np.linspace(0, 10, 4) +mesh.z_grid = np.linspace(0, 5, 5) + +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=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") + + +@pytest.mark.parametrize("this", [True, False]) +def test_rz_slice_of_data_flat_data_normalized(volume_normalization): + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=flatt_data, + axis='R-Z', + slice_index=slice_index, + volume_normalization=volume_normalization, + ) + extent = mesh.get_mpl_plot_extent() + x_label, y_label = mesh.get_axis_labels() + plt.xlabel(x_label) + plt.ylabel(y_label) + plt.imshow(data, extent=extent) + assert data.shape == (4, 3) + + +@pytest.mark.parametrize("this", [True, False]) +def test_rz_slice_of_data_simulation_normalization(volume_normalization): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=my_tally_result.mean, + axis='R-Z', slice_index=slice_index, - volume_normalization=True, + volume_normalization=volume_normalization, ) extent = mesh.get_mpl_plot_extent() x_label, y_label = mesh.get_axis_labels() @@ -70,3 +91,25 @@ def test_slice_of_data(): plt.ylabel(y_label) plt.imshow(data, extent=extent) assert data.shape == (4, 3) + +def test_phiz_slice_of_data_simulation(): + + for slice_index in range(len(mesh.z_grid) - 1): + theta, r, values = mesh.slice_of_data( + dataset=my_tally_result.mean, + slice_index=slice_index, + axis='Phi-R', + volume_normalization=False, + ) + # plt.xlabel(x_label) + # plt.ylabel(y_label) + + fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) + im = ax.contourf(theta, r, values) # , locator=ticker.LogLocator()) + + plt.colorbar(im) + + plt.show() + + + From 5f3aad1b474ac268073b1ef7386e93e959bcb42c Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 23 Mar 2023 12:33:19 +0000 Subject: [PATCH 4/9] [skip ci] Apply formatting changes --- examples/plot_phir_slices.py | 2 +- src/openmc_cylindrical_mesh_plotter/core.py | 10 +++++----- tests/test_slice_of_data.py | 12 ++++-------- 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/examples/plot_phir_slices.py b/examples/plot_phir_slices.py index fdb941e..d5055c3 100644 --- a/examples/plot_phir_slices.py +++ b/examples/plot_phir_slices.py @@ -67,7 +67,7 @@ theta, r, values = mesh.slice_of_data( dataset=my_tally_result.mean, slice_index=slice_index, - axis='Phi-R', + axis="Phi-R", volume_normalization=False, ) # plt.xlabel(x_label) diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index 45fdabb..2817980 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -10,19 +10,19 @@ def slice_of_data( slice_index=0, volume_normalization: bool = True, ): - if axis == 'R-Z': + if axis == "R-Z": return slice_of_rz_data( self, dataset=dataset, slice_index=slice_index, - volume_normalization=volume_normalization + volume_normalization=volume_normalization, ) - elif axis == 'Phi-R': + elif axis == "Phi-R": return slice_of_phir_data( self, dataset=dataset, slice_index=slice_index, - volume_normalization=volume_normalization + volume_normalization=volume_normalization, ) else: raise ValueError(f'axis must be either "R-Z" or "Phi-R", not {axis}') @@ -50,13 +50,13 @@ def slice_of_phir_data( return theta, r, values + def slice_of_rz_data( self, dataset: np.ndarray, slice_index=0, volume_normalization: bool = True, ): - lower_index = int(slice_index * (len(self.r_grid) - 1)) upper_index = int((slice_index + 1) * (len(self.r_grid) - 1)) diff --git a/tests/test_slice_of_data.py b/tests/test_slice_of_data.py index af032ee..47a1b37 100644 --- a/tests/test_slice_of_data.py +++ b/tests/test_slice_of_data.py @@ -63,7 +63,7 @@ def test_rz_slice_of_data_flat_data_normalized(volume_normalization): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flatt_data, - axis='R-Z', + axis="R-Z", slice_index=slice_index, volume_normalization=volume_normalization, ) @@ -77,11 +77,10 @@ def test_rz_slice_of_data_flat_data_normalized(volume_normalization): @pytest.mark.parametrize("this", [True, False]) def test_rz_slice_of_data_simulation_normalization(volume_normalization): - for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=my_tally_result.mean, - axis='R-Z', + axis="R-Z", slice_index=slice_index, volume_normalization=volume_normalization, ) @@ -92,13 +91,13 @@ def test_rz_slice_of_data_simulation_normalization(volume_normalization): plt.imshow(data, extent=extent) assert data.shape == (4, 3) -def test_phiz_slice_of_data_simulation(): +def test_phiz_slice_of_data_simulation(): for slice_index in range(len(mesh.z_grid) - 1): theta, r, values = mesh.slice_of_data( dataset=my_tally_result.mean, slice_index=slice_index, - axis='Phi-R', + axis="Phi-R", volume_normalization=False, ) # plt.xlabel(x_label) @@ -110,6 +109,3 @@ def test_phiz_slice_of_data_simulation(): plt.colorbar(im) plt.show() - - - From 71771156fa01d57e7d5b497b451ee79003230828 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 23 Mar 2023 18:57:27 +0000 Subject: [PATCH 5/9] added tests --- src/openmc_cylindrical_mesh_plotter/core.py | 2 + tests/test_slice_of_data.py | 256 +++++++++++++++----- 2 files changed, 193 insertions(+), 65 deletions(-) diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index 2817980..abf1288 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -42,6 +42,8 @@ def slice_of_phir_data( lower_index = int(slice_index * (len(self.phi_grid) - 1)) upper_index = int((slice_index + 1) * (len(self.phi_grid) - 1)) + + print(dataset.flatten()) # both order A and C appear to work # values=dataset.flatten().reshape(-1,len(self.r_grid)-1,order='C')[:len(self.phi_grid)-1] values = dataset.flatten().reshape(-1, len(self.r_grid) - 1, order="A")[ diff --git a/tests/test_slice_of_data.py b/tests/test_slice_of_data.py index 47a1b37..54d3461 100644 --- a/tests/test_slice_of_data.py +++ b/tests/test_slice_of_data.py @@ -3,6 +3,7 @@ from math import pi import matplotlib.pyplot as plt import openmc_cylindrical_mesh_plotter # adds slice_of_data method to CylindricalMesh +import pytest mesh = openmc.CylindricalMesh() @@ -10,102 +11,227 @@ mesh.r_grid = np.linspace(0, 10, 4) mesh.z_grid = np.linspace(0, 5, 5) -tally = openmc.Tally(name="my_tally") -mesh_filter = openmc.MeshFilter(mesh) -tally.filters.append(mesh_filter) -tally.scores.append("flux") -tallies = openmc.Tallies([tally]) +@pytest.fixture +def circular_source_simulation(): -outer_surface = openmc.Sphere(r=100, boundary_type="vacuum") -cell = openmc.Cell(region=-outer_surface) + tally = openmc.Tally(name="my_tally") + mesh_filter = openmc.MeshFilter(mesh) + tally.filters.append(mesh_filter) + tally.scores.append("flux") + tallies = openmc.Tallies([tally]) -material = openmc.Material() -material.add_nuclide("Fe56", 1) -material.set_density("g/cm3", 0.1) -my_materials = openmc.Materials([material]) + outer_surface = openmc.Sphere(r=100, boundary_type="vacuum") + cell = openmc.Cell(region=-outer_surface) -universe = openmc.Universe(cells=[cell]) -my_geometry = openmc.Geometry(universe) + material = openmc.Material() + material.add_nuclide("Fe56", 1) + material.set_density("g/cm3", 0.1) + my_materials = openmc.Materials([material]) -my_source = openmc.Source() + universe = openmc.Universe(cells=[cell]) + my_geometry = openmc.Geometry(universe) -# 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=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_source = openmc.Source() -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 + # 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=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() -model = openmc.model.Model(my_geometry, my_materials, my_settings, tallies) -sp_filename = model.run() + 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 -statepoint = openmc.StatePoint(sp_filename) + model = openmc.model.Model(my_geometry, my_materials, my_settings, tallies) + sp_filename = model.run() -my_tally_result = statepoint.get_tally(name="my_tally") + statepoint = openmc.StatePoint(sp_filename) + my_tally_result = statepoint.get_tally(name="my_tally") + + return my_tally_result.mean + +@pytest.fixture +def point_source_simulation(): + mesh = openmc.CylindricalMesh() + mesh.phi_grid = np.linspace(0.0, 2 * pi, 10) + mesh.r_grid = np.linspace(0, 10, 4) + mesh.z_grid = np.linspace(0, 5, 5) + + 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() + + my_source.space = openmc.stats.Point((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") + + return my_tally_result.mean + +@pytest.fixture +def flat_data(): + flat_data = [1] * len(mesh.phi_grid) * len(mesh.r_grid) * len(mesh.z_grid) + return np.array(flat_data) + + +def test_get_mpl_plot_extent(): + pass + #todo get extern for both plots polar and imshow + + +def test_get_axis_labels(): + #todo get labels for both plots polar and imshow + pass + +def test_rz_slice_of_data_flat_data_normalized(flat_data): -@pytest.mark.parametrize("this", [True, False]) -def test_rz_slice_of_data_flat_data_normalized(volume_normalization): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( - dataset=flatt_data, + dataset=flat_data, axis="R-Z", slice_index=slice_index, - volume_normalization=volume_normalization, + volume_normalization=True, ) - extent = mesh.get_mpl_plot_extent() - x_label, y_label = mesh.get_axis_labels() - plt.xlabel(x_label) - plt.ylabel(y_label) - plt.imshow(data, extent=extent) + assert data.shape == (4, 3) + #TODO check the data values are smaller on the right than the left + # source is in the middle, voxel volumes are smaller in the center + +def test_rz_slice_of_data_flat_data_unnormalized(flat_data): -@pytest.mark.parametrize("this", [True, False]) -def test_rz_slice_of_data_simulation_normalization(volume_normalization): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( - dataset=my_tally_result.mean, + dataset=flat_data, axis="R-Z", slice_index=slice_index, - volume_normalization=volume_normalization, + volume_normalization=False, ) - extent = mesh.get_mpl_plot_extent() - x_label, y_label = mesh.get_axis_labels() - plt.xlabel(x_label) - plt.ylabel(y_label) - plt.imshow(data, extent=extent) + assert data.shape == (4, 3) + #TODO check the data values equal everywhere -def test_phiz_slice_of_data_simulation(): - for slice_index in range(len(mesh.z_grid) - 1): - theta, r, values = mesh.slice_of_data( - dataset=my_tally_result.mean, +def test_phir_slice_of_data_flat_data_normalized(flat_data): + + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=flat_data, + axis="Phi-R", + slice_index=slice_index, + volume_normalization=True, + ) + + assert data.shape == (4, 3) + + #TODO check the data values + +def test_phir_slice_of_data_flat_data_unnormalized(flat_data): + + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=flat_data, + axis="R-Z", slice_index=slice_index, + volume_normalization=False, + ) + + assert data.shape == (4, 3) + + #TODO check the data values equal everywhere + +def test_rz_slice_of_data_point_simulation_normalization(point_source_simulation): + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=point_source_simulation, + axis="R-Z", + slice_index=slice_index, + volume_normalization=True, + ) + + assert data.shape == (4, 3) + + #TODO test + +def test_phir_slice_of_data_circular_simulation_normalization(circular_source_simulation): + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=circular_source_simulation, axis="Phi-R", + slice_index=slice_index, + volume_normalization=True, + ) + + assert data.shape == (4, 3) + + #TODO test + +def test_rz_slice_of_data_point_simulation_unnormalization(point_source_simulation): + for slice_index in range(len(mesh.phi_grid) - 1): + data = mesh.slice_of_data( + dataset=point_source_simulation, + axis="R-Z", + slice_index=slice_index, volume_normalization=False, ) - # plt.xlabel(x_label) - # plt.ylabel(y_label) - fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) - im = ax.contourf(theta, r, values) # , locator=ticker.LogLocator()) + assert data.shape == (4, 3) + + #TODO test - plt.colorbar(im) +def test_phir_slice_of_data_circular_simulation_unnormalization(circular_source_simulation): + for slice_index in range(len(mesh.phi_grid) - 1): + theta, r, values = mesh.slice_of_data( + dataset=circular_source_simulation, + axis="Phi-R", + slice_index=slice_index, + volume_normalization=False, + ) + theta.shape == (4, 3) + r.shape == (4, 3) + values.shape == (4, 3) - plt.show() + #TODO test From 19f9ecad4f7b675695a94399e11ec42b605f3ccd Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 23 Mar 2023 18:59:02 +0000 Subject: [PATCH 6/9] [skip ci] Apply formatting changes --- src/openmc_cylindrical_mesh_plotter/core.py | 1 - tests/test_slice_of_data.py | 47 ++++++++++++--------- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index abf1288..f20aeff 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -42,7 +42,6 @@ def slice_of_phir_data( lower_index = int(slice_index * (len(self.phi_grid) - 1)) upper_index = int((slice_index + 1) * (len(self.phi_grid) - 1)) - print(dataset.flatten()) # both order A and C appear to work # values=dataset.flatten().reshape(-1,len(self.r_grid)-1,order='C')[:len(self.phi_grid)-1] diff --git a/tests/test_slice_of_data.py b/tests/test_slice_of_data.py index 54d3461..26619f9 100644 --- a/tests/test_slice_of_data.py +++ b/tests/test_slice_of_data.py @@ -11,9 +11,9 @@ mesh.r_grid = np.linspace(0, 10, 4) mesh.z_grid = np.linspace(0, 5, 5) + @pytest.fixture def circular_source_simulation(): - tally = openmc.Tally(name="my_tally") mesh_filter = openmc.MeshFilter(mesh) tally.filters.append(mesh_filter) @@ -31,7 +31,6 @@ def circular_source_simulation(): universe = openmc.Universe(cells=[cell]) my_geometry = openmc.Geometry(universe) - my_source = openmc.Source() # the distribution of radius is just a single value @@ -64,6 +63,7 @@ def circular_source_simulation(): return my_tally_result.mean + @pytest.fixture def point_source_simulation(): mesh = openmc.CylindricalMesh() @@ -90,7 +90,7 @@ def point_source_simulation(): my_source = openmc.Source() - my_source.space = openmc.stats.Point((0,0., 0)) + my_source.space = openmc.stats.Point((0, 0.0, 0)) # sets the direction to isotropic my_source.angle = openmc.stats.Isotropic() @@ -111,6 +111,7 @@ def point_source_simulation(): return my_tally_result.mean + @pytest.fixture def flat_data(): flat_data = [1] * len(mesh.phi_grid) * len(mesh.r_grid) * len(mesh.z_grid) @@ -119,15 +120,15 @@ def flat_data(): def test_get_mpl_plot_extent(): pass - #todo get extern for both plots polar and imshow + # todo get extern for both plots polar and imshow def test_get_axis_labels(): - #todo get labels for both plots polar and imshow + # todo get labels for both plots polar and imshow pass -def test_rz_slice_of_data_flat_data_normalized(flat_data): +def test_rz_slice_of_data_flat_data_normalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, @@ -138,11 +139,11 @@ def test_rz_slice_of_data_flat_data_normalized(flat_data): assert data.shape == (4, 3) - #TODO check the data values are smaller on the right than the left + # TODO check the data values are smaller on the right than the left # source is in the middle, voxel volumes are smaller in the center -def test_rz_slice_of_data_flat_data_unnormalized(flat_data): +def test_rz_slice_of_data_flat_data_unnormalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, @@ -153,10 +154,10 @@ def test_rz_slice_of_data_flat_data_unnormalized(flat_data): assert data.shape == (4, 3) - #TODO check the data values equal everywhere + # TODO check the data values equal everywhere -def test_phir_slice_of_data_flat_data_normalized(flat_data): +def test_phir_slice_of_data_flat_data_normalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, @@ -167,10 +168,10 @@ def test_phir_slice_of_data_flat_data_normalized(flat_data): assert data.shape == (4, 3) - #TODO check the data values + # TODO check the data values -def test_phir_slice_of_data_flat_data_unnormalized(flat_data): +def test_phir_slice_of_data_flat_data_unnormalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, @@ -181,7 +182,8 @@ def test_phir_slice_of_data_flat_data_unnormalized(flat_data): assert data.shape == (4, 3) - #TODO check the data values equal everywhere + # TODO check the data values equal everywhere + def test_rz_slice_of_data_point_simulation_normalization(point_source_simulation): for slice_index in range(len(mesh.phi_grid) - 1): @@ -194,9 +196,12 @@ def test_rz_slice_of_data_point_simulation_normalization(point_source_simulation assert data.shape == (4, 3) - #TODO test + # TODO test + -def test_phir_slice_of_data_circular_simulation_normalization(circular_source_simulation): +def test_phir_slice_of_data_circular_simulation_normalization( + circular_source_simulation, +): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=circular_source_simulation, @@ -207,7 +212,8 @@ def test_phir_slice_of_data_circular_simulation_normalization(circular_source_si assert data.shape == (4, 3) - #TODO test + # TODO test + def test_rz_slice_of_data_point_simulation_unnormalization(point_source_simulation): for slice_index in range(len(mesh.phi_grid) - 1): @@ -220,9 +226,12 @@ def test_rz_slice_of_data_point_simulation_unnormalization(point_source_simulati assert data.shape == (4, 3) - #TODO test + # TODO test + -def test_phir_slice_of_data_circular_simulation_unnormalization(circular_source_simulation): +def test_phir_slice_of_data_circular_simulation_unnormalization( + circular_source_simulation, +): for slice_index in range(len(mesh.phi_grid) - 1): theta, r, values = mesh.slice_of_data( dataset=circular_source_simulation, @@ -234,4 +243,4 @@ def test_phir_slice_of_data_circular_simulation_unnormalization(circular_source_ r.shape == (4, 3) values.shape == (4, 3) - #TODO test + # TODO test From eb6d2628a25e8b98f71a90ac822cfd248b03fa9c Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 24 Mar 2023 23:17:24 +0000 Subject: [PATCH 7/9] moved to pandas method --- .gitignore | 1 + examples/plot_phir_slice_point_source.py | 72 ++++++++++++ ...ices.py => plot_phir_slice_ring_source.py} | 20 ++-- ...ices.py => plot_rz_slices_point_source.py} | 9 +- examples/plot_rz_slices_ring_source.py | 79 +++++++++++++ src/openmc_cylindrical_mesh_plotter/core.py | 108 +++++++++++++----- tests/test_slice_of_data.py | 20 ++-- 7 files changed, 260 insertions(+), 49 deletions(-) create mode 100644 examples/plot_phir_slice_point_source.py rename examples/{plot_phir_slices.py => plot_phir_slice_ring_source.py} (80%) rename examples/{plot_rz_slices.py => plot_rz_slices_point_source.py} (92%) create mode 100644 examples/plot_rz_slices_ring_source.py diff --git a/.gitignore b/.gitignore index 58d6fc5..8f18586 100644 --- a/.gitignore +++ b/.gitignore @@ -132,3 +132,4 @@ src/*/_version.py *.vtk *.h5 *.out +*.png diff --git a/examples/plot_phir_slice_point_source.py b/examples/plot_phir_slice_point_source.py new file mode 100644 index 0000000..50bc322 --- /dev/null +++ b/examples/plot_phir_slice_point_source.py @@ -0,0 +1,72 @@ +# 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() diff --git a/examples/plot_phir_slices.py b/examples/plot_phir_slice_ring_source.py similarity index 80% rename from examples/plot_phir_slices.py rename to examples/plot_phir_slice_ring_source.py index d5055c3..b957e40 100644 --- a/examples/plot_phir_slices.py +++ b/examples/plot_phir_slice_ring_source.py @@ -8,9 +8,10 @@ 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, 2 * pi, 10) +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) @@ -38,7 +39,7 @@ # 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=1 * 3.14159265359) +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)) @@ -63,19 +64,20 @@ my_tally_result = statepoint.get_tally(name="my_tally") -for slice_index in range(len(mesh.z_grid) - 1): +for slice_index in range(1, len(mesh.z_grid)): theta, r, values = mesh.slice_of_data( - dataset=my_tally_result.mean, + dataset=my_tally_result.mean.flatten(), slice_index=slice_index, - axis="Phi-R", + axis="PhiR", volume_normalization=False, ) - # plt.xlabel(x_label) - # plt.ylabel(y_label) fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) - im = ax.contourf(theta, r, values) # , locator=ticker.LogLocator()) + 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) + plt.colorbar(im, label='Flux') plt.show() diff --git a/examples/plot_rz_slices.py b/examples/plot_rz_slices_point_source.py similarity index 92% rename from examples/plot_rz_slices.py rename to examples/plot_rz_slices_point_source.py index db7fb0e..587f28d 100644 --- a/examples/plot_rz_slices.py +++ b/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() diff --git a/examples/plot_rz_slices_ring_source.py b/examples/plot_rz_slices_ring_source.py new file mode 100644 index 0000000..b8d104a --- /dev/null +++ b/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() diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index f20aeff..1727914 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -1,5 +1,6 @@ import matplotlib.pyplot as plt import numpy as np +import pandas as pd import openmc @@ -10,14 +11,14 @@ def slice_of_data( slice_index=0, volume_normalization: bool = True, ): - if axis == "R-Z": + if axis == "RZ": return slice_of_rz_data( self, dataset=dataset, slice_index=slice_index, volume_normalization=volume_normalization, ) - elif axis == "Phi-R": + elif axis == "PhiR": return slice_of_phir_data( self, dataset=dataset, @@ -25,7 +26,7 @@ def slice_of_data( volume_normalization=volume_normalization, ) else: - raise ValueError(f'axis must be either "R-Z" or "Phi-R", not {axis}') + raise ValueError(f'axis must be either "RZ" or "PhiR", not {axis}') def slice_of_phir_data( @@ -37,19 +38,45 @@ def slice_of_phir_data( actual = np.linspace(self.phi_grid[0], self.phi_grid[-1], self.dimension[1]) expected = np.linspace(self.r_grid[0], self.r_grid[-1], self.dimension[0]) - r, theta = np.meshgrid(expected, actual) + rad, theta = np.meshgrid(expected, actual) + + # old method which suffered from errors reshaping but avoided using pandas + # lower_index = int(slice_index * (len(self.phi_grid) - 1)) + # upper_index = int((slice_index + 1) * (len(self.phi_grid) - 1)) + + # print(dataset.flatten()) + # # both order A and C appear to work + # # values=dataset.flatten().reshape(-1,len(self.r_grid)-1,order='C')[:len(self.phi_grid)-1] + # values = dataset.flatten().reshape(-1, len(self.r_grid) - 1, order="A")[ + # lower_index:upper_index + # ] + + # return theta, r, values + indices = list(self.indices) + + r = [entry[0] for entry in indices] + phi = [entry[1] for entry in indices] + z = [entry[2] for entry in indices] + + df = pd.DataFrame( + { + 'r': r, + 'phi': phi, + 'z': z, + 'value': dataset, + 'volume':self.volumes.flatten(), + 'normalised_value':dataset/self.volumes.flatten(), + } + ) + df_slice = df[df['z'] == slice_index] - lower_index = int(slice_index * (len(self.phi_grid) - 1)) - upper_index = int((slice_index + 1) * (len(self.phi_grid) - 1)) - - print(dataset.flatten()) - # both order A and C appear to work - # values=dataset.flatten().reshape(-1,len(self.r_grid)-1,order='C')[:len(self.phi_grid)-1] - values = dataset.flatten().reshape(-1, len(self.r_grid) - 1, order="A")[ - lower_index:upper_index - ] + if volume_normalization: + shaped_slice = df_slice['normalised_value'].to_numpy().reshape(-1, len(self.r_grid)-1) + else: + shaped_slice = df_slice['value'].to_numpy().reshape(-1, len(self.r_grid)-1) + + return theta, rad, shaped_slice - return theta, r, values def slice_of_rz_data( @@ -58,22 +85,49 @@ def slice_of_rz_data( slice_index=0, volume_normalization: bool = True, ): - lower_index = int(slice_index * (len(self.r_grid) - 1)) - upper_index = int((slice_index + 1) * (len(self.r_grid) - 1)) + # old method which suffered from errors reshaping but avoided using pandas + # lower_index = int(slice_index * (len(self.r_grid) - 1)) + # upper_index = int((slice_index + 1) * (len(self.r_grid) - 1)) + + # if volume_normalization: + # data_slice = dataset.flatten().reshape(-1, len(self.z_grid) - 1, order="F")[ + # lower_index:upper_index + # ] + # data_slice = data_slice / self.volumes[:, 1, :] + + # return np.rot90(data_slice) + + # data_slice = dataset.flatten().reshape(-1, len(self.z_grid) - 1, order="F")[ + # lower_index:upper_index + # ] + + # return np.rot90(data_slice) + + indices = list(self.indices) + + r = [entry[0] for entry in indices] + phi = [entry[1] for entry in indices] + z = [entry[2] for entry in indices] + + df = pd.DataFrame( + { + 'r': r, + 'phi': phi, + 'z': z, + 'value': dataset, + 'volume':self.volumes.flatten(), + 'normalised_value':dataset/self.volumes.flatten(), + } + ) + df_slice = df[df['phi'] == slice_index] if volume_normalization: - data_slice = dataset.flatten().reshape(-1, len(self.z_grid) - 1, order="F")[ - lower_index:upper_index - ] - data_slice = data_slice / self.volumes[:, 1, :] - - return np.rot90(data_slice) - - data_slice = dataset.flatten().reshape(-1, len(self.z_grid) - 1, order="F")[ - lower_index:upper_index - ] + shaped_slice = df_slice['normalised_value'].to_numpy().reshape(-1, len(self.r_grid)-1) + else: + shaped_slice = df_slice['value'].to_numpy().reshape(-1, len(self.r_grid)-1) + + return np.flipud(shaped_slice) - return np.rot90(data_slice) def get_mpl_plot_extent(self): diff --git a/tests/test_slice_of_data.py b/tests/test_slice_of_data.py index 26619f9..ee31dbc 100644 --- a/tests/test_slice_of_data.py +++ b/tests/test_slice_of_data.py @@ -61,7 +61,7 @@ def circular_source_simulation(): my_tally_result = statepoint.get_tally(name="my_tally") - return my_tally_result.mean + return my_tally_result.mean.flatten() @pytest.fixture @@ -109,7 +109,7 @@ def point_source_simulation(): my_tally_result = statepoint.get_tally(name="my_tally") - return my_tally_result.mean + return my_tally_result.mean.flatten() @pytest.fixture @@ -132,7 +132,7 @@ def test_rz_slice_of_data_flat_data_normalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, - axis="R-Z", + axis="RZ", slice_index=slice_index, volume_normalization=True, ) @@ -147,7 +147,7 @@ def test_rz_slice_of_data_flat_data_unnormalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, - axis="R-Z", + axis="RZ", slice_index=slice_index, volume_normalization=False, ) @@ -161,7 +161,7 @@ def test_phir_slice_of_data_flat_data_normalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, - axis="Phi-R", + axis="PhiR", slice_index=slice_index, volume_normalization=True, ) @@ -175,7 +175,7 @@ def test_phir_slice_of_data_flat_data_unnormalized(flat_data): for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=flat_data, - axis="R-Z", + axis="RZ", slice_index=slice_index, volume_normalization=False, ) @@ -189,7 +189,7 @@ def test_rz_slice_of_data_point_simulation_normalization(point_source_simulation for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=point_source_simulation, - axis="R-Z", + axis="RZ", slice_index=slice_index, volume_normalization=True, ) @@ -205,7 +205,7 @@ def test_phir_slice_of_data_circular_simulation_normalization( for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=circular_source_simulation, - axis="Phi-R", + axis="PhiR", slice_index=slice_index, volume_normalization=True, ) @@ -219,7 +219,7 @@ def test_rz_slice_of_data_point_simulation_unnormalization(point_source_simulati for slice_index in range(len(mesh.phi_grid) - 1): data = mesh.slice_of_data( dataset=point_source_simulation, - axis="R-Z", + axis="RZ", slice_index=slice_index, volume_normalization=False, ) @@ -235,7 +235,7 @@ def test_phir_slice_of_data_circular_simulation_unnormalization( for slice_index in range(len(mesh.phi_grid) - 1): theta, r, values = mesh.slice_of_data( dataset=circular_source_simulation, - axis="Phi-R", + axis="PhiR", slice_index=slice_index, volume_normalization=False, ) From 2184556d020df8c239f2a758f83ef1e455436a72 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 24 Mar 2023 23:17:48 +0000 Subject: [PATCH 8/9] [skip ci] Apply formatting changes --- examples/plot_phir_slice_point_source.py | 8 ++-- examples/plot_phir_slice_ring_source.py | 10 +++-- examples/plot_rz_slices_point_source.py | 8 ++-- examples/plot_rz_slices_ring_source.py | 8 ++-- src/openmc_cylindrical_mesh_plotter/core.py | 46 +++++++++++---------- 5 files changed, 44 insertions(+), 36 deletions(-) diff --git a/examples/plot_phir_slice_point_source.py b/examples/plot_phir_slice_point_source.py index 50bc322..3efc71a 100644 --- a/examples/plot_phir_slice_point_source.py +++ b/examples/plot_phir_slice_point_source.py @@ -34,7 +34,7 @@ my_source = openmc.Source() # this makes a point source instead with -my_source.space = openmc.stats.Point((0,-5,0)) +my_source.space = openmc.stats.Point((0, -5, 0)) # sets the direction to isotropic my_source.angle = openmc.stats.Isotropic() @@ -62,11 +62,13 @@ ) fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) - im = ax.contourf(theta, r, values, extent=(0,100,0,50)) # , locator=ticker.LogLocator()) + 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.colorbar(im, label="Flux") plt.show() diff --git a/examples/plot_phir_slice_ring_source.py b/examples/plot_phir_slice_ring_source.py index b957e40..32c63bc 100644 --- a/examples/plot_phir_slice_ring_source.py +++ b/examples/plot_phir_slice_ring_source.py @@ -11,7 +11,9 @@ 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.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) @@ -73,11 +75,13 @@ ) fig, ax = plt.subplots(subplot_kw=dict(projection="polar")) - im = ax.contourf(theta, r, values, extent=(0,100,0,50)) # , locator=ticker.LogLocator()) + 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.colorbar(im, label="Flux") plt.show() diff --git a/examples/plot_rz_slices_point_source.py b/examples/plot_rz_slices_point_source.py index 587f28d..b97b4ec 100644 --- a/examples/plot_rz_slices_point_source.py +++ b/examples/plot_rz_slices_point_source.py @@ -61,10 +61,10 @@ my_tally_result = statepoint.get_tally(name="my_tally") -for slice_index in range(1,len(mesh.phi_grid)): +for slice_index in range(1, len(mesh.phi_grid)): data = mesh.slice_of_data( dataset=my_tally_result.mean.flatten(), - axis='RZ', + axis="RZ", # dataset=np.array(2*19*49*[1]), flat data for testing slice_index=slice_index, volume_normalization=False, @@ -74,7 +74,7 @@ plt.xlabel(x_label) plt.ylabel(y_label) im = plt.imshow(data, extent=extent) - - plt.colorbar(im, label='Flux') + + plt.colorbar(im, label="Flux") plt.show() diff --git a/examples/plot_rz_slices_ring_source.py b/examples/plot_rz_slices_ring_source.py index b8d104a..654ae06 100644 --- a/examples/plot_rz_slices_ring_source.py +++ b/examples/plot_rz_slices_ring_source.py @@ -61,19 +61,19 @@ my_tally_result = statepoint.get_tally(name="my_tally") -for slice_index in range(1,len(mesh.phi_grid)): +for slice_index in range(1, len(mesh.phi_grid)): data = mesh.slice_of_data( dataset=my_tally_result.mean.flatten(), - axis='RZ', + 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.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.colorbar(im, label="Flux") plt.show() diff --git a/src/openmc_cylindrical_mesh_plotter/core.py b/src/openmc_cylindrical_mesh_plotter/core.py index 1727914..c4d515d 100644 --- a/src/openmc_cylindrical_mesh_plotter/core.py +++ b/src/openmc_cylindrical_mesh_plotter/core.py @@ -60,23 +60,24 @@ def slice_of_phir_data( df = pd.DataFrame( { - 'r': r, - 'phi': phi, - 'z': z, - 'value': dataset, - 'volume':self.volumes.flatten(), - 'normalised_value':dataset/self.volumes.flatten(), + "r": r, + "phi": phi, + "z": z, + "value": dataset, + "volume": self.volumes.flatten(), + "normalised_value": dataset / self.volumes.flatten(), } ) - df_slice = df[df['z'] == slice_index] + df_slice = df[df["z"] == slice_index] if volume_normalization: - shaped_slice = df_slice['normalised_value'].to_numpy().reshape(-1, len(self.r_grid)-1) + shaped_slice = ( + df_slice["normalised_value"].to_numpy().reshape(-1, len(self.r_grid) - 1) + ) else: - shaped_slice = df_slice['value'].to_numpy().reshape(-1, len(self.r_grid)-1) - - return theta, rad, shaped_slice + shaped_slice = df_slice["value"].to_numpy().reshape(-1, len(self.r_grid) - 1) + return theta, rad, shaped_slice def slice_of_rz_data( @@ -111,23 +112,24 @@ def slice_of_rz_data( df = pd.DataFrame( { - 'r': r, - 'phi': phi, - 'z': z, - 'value': dataset, - 'volume':self.volumes.flatten(), - 'normalised_value':dataset/self.volumes.flatten(), + "r": r, + "phi": phi, + "z": z, + "value": dataset, + "volume": self.volumes.flatten(), + "normalised_value": dataset / self.volumes.flatten(), } ) - df_slice = df[df['phi'] == slice_index] + df_slice = df[df["phi"] == slice_index] if volume_normalization: - shaped_slice = df_slice['normalised_value'].to_numpy().reshape(-1, len(self.r_grid)-1) + shaped_slice = ( + df_slice["normalised_value"].to_numpy().reshape(-1, len(self.r_grid) - 1) + ) else: - shaped_slice = df_slice['value'].to_numpy().reshape(-1, len(self.r_grid)-1) - - return np.flipud(shaped_slice) + shaped_slice = df_slice["value"].to_numpy().reshape(-1, len(self.r_grid) - 1) + return np.flipud(shaped_slice) def get_mpl_plot_extent(self): From a77a87f64275ccbd669c4d40771c62eeaca0c499 Mon Sep 17 00:00:00 2001 From: shimwell Date: Fri, 24 Mar 2023 23:21:10 +0000 Subject: [PATCH 9/9] added images to readme --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index cf8b84f..6de7911 100644 --- a/README.md +++ b/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