From facc1995e342fbfd8e0c13ff42550beabec5d4fe Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 5 Oct 2022 09:29:19 +0200 Subject: [PATCH 01/58] marks --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index a4c50a6084..4059fee2d3 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -165,7 +165,7 @@ class StructuredMesh(MeshBase): pass @property - def vertices(self): + def vertices(self): # TODO needs changing """Return coordinates of mesh vertices. Returns @@ -178,7 +178,7 @@ class StructuredMesh(MeshBase): return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=-1) @property - def centroids(self): + def centroids(self): # TODO needs changing """Return coordinates of mesh element centroids. Returns From 2a4c8bbb9b0253562131ddd207ebf1be523edc69 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 06:47:11 +0000 Subject: [PATCH 02/58] centre attribute to CylindricalMesh and SphericalMesh --- openmc/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 4059fee2d3..59d64a1464 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,6 +1052,7 @@ class CylindricalMesh(StructuredMesh): self._r_grid = None self._phi_grid = [0.0, 2*pi] self._z_grid = None + self._centre = [0, 0, 0] @property def dimension(self): @@ -1284,6 +1285,7 @@ class SphericalMesh(StructuredMesh): self._r_grid = None self._theta_grid = [0, pi] self._phi_grid = [0, 2*pi] + self._centre = [0, 0, 0] @property def dimension(self): From c20b351a2fa26d8198f833b256d140178b4c4e9c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 06:49:13 +0000 Subject: [PATCH 03/58] write to xml --- openmc/mesh.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 59d64a1464..cc828acdc0 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1161,6 +1161,9 @@ class CylindricalMesh(StructuredMesh): subelement = ET.SubElement(element, "z_grid") subelement.text = ' '.join(map(str, self.z_grid)) + subelement = ET.SubElement(element, "centre") + subelement.text = ' '.join(map(str, self.centre)) + return element @classmethod @@ -1184,6 +1187,8 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] + # TODO: add read centre + return mesh @property @@ -1394,6 +1399,9 @@ class SphericalMesh(StructuredMesh): subelement = ET.SubElement(element, "phi_grid") subelement.text = ' '.join(map(str, self.phi_grid)) + subelement = ET.SubElement(element, "centre") + subelement.text = ' '.join(map(str, self.centre)) + return element @classmethod @@ -1417,6 +1425,8 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] + + # TODO: add read centre return mesh @property From 81ac9b0fd900df2d21cb194fbf242375844c9658 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 06:50:45 +0000 Subject: [PATCH 04/58] setters/getters --- openmc/mesh.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index cc828acdc0..3efa1e1149 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1064,6 +1064,10 @@ class CylindricalMesh(StructuredMesh): def n_dimension(self): return 3 + @property + def centre(self): + return self._centre + @property def r_grid(self): return self._r_grid @@ -1090,6 +1094,11 @@ class CylindricalMesh(StructuredMesh): for p in range(1, np + 1) for r in range(1, nr + 1)) + @centre.setter + def centre(self, coords): + cv.check_type('mesh r_grid', coords, Iterable, Real) + self._centre = np.asarray(coords) + @r_grid.setter def r_grid(self, grid): cv.check_type('mesh r_grid', grid, Iterable, Real) @@ -1302,6 +1311,10 @@ class SphericalMesh(StructuredMesh): def n_dimension(self): return 3 + @property + def centre(self): + return self._centre + @property def r_grid(self): return self._r_grid @@ -1328,6 +1341,11 @@ class SphericalMesh(StructuredMesh): for t in range(1, nt + 1) for r in range(1, nr + 1)) + @centre.setter + def centre(self, coords): + cv.check_type('mesh r_grid', coords, Iterable, Real) + self._centre = np.asarray(coords) + @r_grid.setter def r_grid(self, grid): cv.check_type('mesh r_grid', grid, Iterable, Real) From f0f5e34b33e79fe5818ae227bd5da192a208bf40 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 06:55:56 +0000 Subject: [PATCH 05/58] read centre from xml --- src/mesh.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index d9ad0ce8ad..8307e5fe87 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,6 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); + centre = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -983,6 +984,11 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( mapped_r[1] += 2 * M_PI; } + // TODO: pass centre as argument + mapped_r[0] += centre[0] + mapped_r[1] += centre[1] + mapped_r[2] += centre[2] + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_phi(idx[1]); @@ -1187,6 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); + centre = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1216,6 +1223,11 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[2] += 2 * M_PI; } + // TODO: pass centre as argument + mapped_r[0] += centre[0] + mapped_r[1] += centre[1] + mapped_r[2] += centre[2] + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_theta(idx[1]); From 5787616ae41593bac4d100da162da5a8f98a1cea Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 07:18:00 +0000 Subject: [PATCH 06/58] added two tests to check centre was written/read to/from xml --- tests/unit_tests/test_cylindrical_mesh.py | 30 +++++++++++++++++++++++ tests/unit_tests/test_spherical_mesh.py | 30 +++++++++++++++++++++++ 2 files changed, 60 insertions(+) create mode 100644 tests/unit_tests/test_cylindrical_mesh.py create mode 100644 tests/unit_tests/test_spherical_mesh.py diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py new file mode 100644 index 0000000000..89322d7045 --- /dev/null +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -0,0 +1,30 @@ +import openmc +import numpy as np + +def test_centre_read_write_to_xml(): + """Tests that the centre attribute can be written and read back to XML + """ + # build + mesh = openmc.CylindricalMesh() + mesh.phi_grid = [1, 2, 3] + mesh.z_grid = [1, 2, 3] + mesh.r_grid = [1, 2, 3] + + mesh.centre = [0.1, 0.2, 0.3] + + tally = openmc.Tally() + + mesh_filter = openmc.MeshFilter(mesh) + tally.filters.append(mesh_filter) + + tally.scores.append("heating") + + tallies = openmc.Tallies([tally]) + + tallies.export_to_xml() + + # read back + new_tallies = openmc.Tallies.from_xml() + new_tally = new_tallies[0] + new_mesh = new_tally.filters[0].mesh + assert np.allclose(new_mesh.centre, mesh.centre) \ No newline at end of file diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py new file mode 100644 index 0000000000..cdeb188ece --- /dev/null +++ b/tests/unit_tests/test_spherical_mesh.py @@ -0,0 +1,30 @@ +import openmc +import numpy as np + +def test_centre_read_write_to_xml(): + """Tests that the centre attribute can be written and read back to XML + """ + # build + mesh = openmc.SphericalMesh() + mesh.phi_grid = [1, 2, 3] + mesh.theta_grid = [1, 2, 3] + mesh.r_grid = [1, 2, 3] + + mesh.centre = [0.1, 0.2, 0.3] + + tally = openmc.Tally() + + mesh_filter = openmc.MeshFilter(mesh) + tally.filters.append(mesh_filter) + + tally.scores.append("heating") + + tallies = openmc.Tallies([tally]) + + tallies.export_to_xml() + + # read back + new_tallies = openmc.Tallies.from_xml() + new_tally = new_tallies[0] + new_mesh = new_tally.filters[0].mesh + assert np.allclose(new_mesh.centre, mesh.centre) \ No newline at end of file From 16add87c05a2a228daabb4f99a3c09c0defe2b46 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 07:18:52 +0000 Subject: [PATCH 07/58] now reads back from xml --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3efa1e1149..abc46c70b4 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1196,7 +1196,7 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] - # TODO: add read centre + mesh.centre = [float(x) for x in get_text(elem, "centre").split()] return mesh @@ -1443,8 +1443,8 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] + mesh.centre = [float(x) for x in get_text(elem, "centre").split()] - # TODO: add read centre return mesh @property From 21200d9b6eb6f07c0dbb1e2885c6c367dea54100 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 07:42:34 +0000 Subject: [PATCH 08/58] offsets the write_to_vtk coordinates --- openmc/mesh.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index abc46c70b4..8f47e4e7b1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,7 +1052,7 @@ class CylindricalMesh(StructuredMesh): self._r_grid = None self._phi_grid = [0.0, 2*pi] self._z_grid = None - self._centre = [0, 0, 0] + self._centre = [0., 0., 0.] @property def dimension(self): @@ -1251,6 +1251,11 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) + # offset with centre + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2] + return super().write_data_to_vtk( points=pts_cartesian, filename=filename, @@ -1299,7 +1304,7 @@ class SphericalMesh(StructuredMesh): self._r_grid = None self._theta_grid = [0, pi] self._phi_grid = [0, 2*pi] - self._centre = [0, 0, 0] + self._centre = [0., 0., 0.] @property def dimension(self): @@ -1500,6 +1505,11 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) pts_cartesian[:, 2] = r * np.cos(phi) + # offset with centre + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2] + return super().write_data_to_vtk( points=pts_cartesian, filename=filename, From 896ae2792ec6e1dfdaa2f7054dda5cb890d7e43f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 08:40:04 +0000 Subject: [PATCH 09/58] centre to centre_ --- src/mesh.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 8307e5fe87..5fa46cd855 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - centre = get_node_array(node, "centre"); + centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -985,9 +985,9 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( } // TODO: pass centre as argument - mapped_r[0] += centre[0] - mapped_r[1] += centre[1] - mapped_r[2] += centre[2] + mapped_r[0] += centre_[0] + mapped_r[1] += centre_[1] + mapped_r[2] += centre_[2] MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - centre = get_node_array(node, "centre"); + centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1224,9 +1224,9 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( } // TODO: pass centre as argument - mapped_r[0] += centre[0] - mapped_r[1] += centre[1] - mapped_r[2] += centre[2] + mapped_r[0] += centre_[0] + mapped_r[1] += centre_[1] + mapped_r[2] += centre_[2] MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); From d194e000803acc3c6e35c9e95d27859ae43d6b99 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 09:00:27 +0000 Subject: [PATCH 10/58] test if grid_ is defined in this scope --- src/mesh.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index 5fa46cd855..7a0a243f04 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1223,6 +1223,9 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[2] += 2 * M_PI; } + // test + coucou = grid_[0] + // TODO: pass centre as argument mapped_r[0] += centre_[0] mapped_r[1] += centre_[1] From a104c1f890055f5822d2b93db08bdcad255a17a0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 09:07:51 +0000 Subject: [PATCH 11/58] [skip ci] removed test --- src/mesh.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 7a0a243f04..5fa46cd855 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1223,9 +1223,6 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[2] += 2 * M_PI; } - // test - coucou = grid_[0] - // TODO: pass centre as argument mapped_r[0] += centre_[0] mapped_r[1] += centre_[1] From bdd13ba3f4c29512f0fc24d4cf856f010a982e1c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 09:15:27 +0000 Subject: [PATCH 12/58] attempt #2 --- src/mesh.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 5fa46cd855..9fde3dbe15 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - centre_ = get_node_array(node, "centre"); + double centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - centre_ = get_node_array(node, "centre"); + double centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); From e394637f2f873ca3cec6417a02031f340650e368 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 09:32:08 +0000 Subject: [PATCH 13/58] declare variable with auto --- src/mesh.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 9fde3dbe15..ad67b2ba41 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - double centre_ = get_node_array(node, "centre"); + auto centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - double centre_ = get_node_array(node, "centre"); + auto centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); From 61b4f30dd4f1ce0bd44e339b59a7cfe31fef057f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 10:59:35 +0000 Subject: [PATCH 14/58] declaring centre_ in C file --- include/openmc/mesh.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index eece0f1ff1..e6f5ce98b9 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,6 +359,8 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; + array centre_; + int set_grid(); @@ -413,6 +415,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; + array centre_; int set_grid(); From fc5aa68d732fb15e99e103d1543451dbb57f99ab Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 10:59:52 +0000 Subject: [PATCH 15/58] removed auto + ; --- src/mesh.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index ad67b2ba41..c049549eb2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - auto centre_ = get_node_array(node, "centre"); + centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -985,9 +985,9 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( } // TODO: pass centre as argument - mapped_r[0] += centre_[0] - mapped_r[1] += centre_[1] - mapped_r[2] += centre_[2] + mapped_r[0] += centre_[0]; + mapped_r[1] += centre_[1]; + mapped_r[2] += centre_[2]; MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - auto centre_ = get_node_array(node, "centre"); + centre_ = get_node_array(node, "centre"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1224,9 +1224,9 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( } // TODO: pass centre as argument - mapped_r[0] += centre_[0] - mapped_r[1] += centre_[1] - mapped_r[2] += centre_[2] + mapped_r[0] += centre_[0]; + mapped_r[1] += centre_[1]; + mapped_r[2] += centre_[2]; MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); From e1dbe49bbbcae3c227190f48e2d6d89945b64187 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 11:06:24 +0000 Subject: [PATCH 16/58] correct type for centre_ --- include/openmc/mesh.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index e6f5ce98b9..ea6b01aa60 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,7 +359,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array centre_; + vector centre_; int set_grid(); @@ -415,7 +415,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array centre_; + vector centre_; int set_grid(); From 1ed096f6608b0caad5f2f17f4cfc34c278d78b95 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 7 Oct 2022 12:53:34 +0000 Subject: [PATCH 17/58] added centre tag to standard --- tests/regression_tests/filter_mesh/inputs_true.dat | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index d0c481d604..c42f3cd439 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -57,11 +57,13 @@ 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586 -7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5 + 0.0 0.0 0.0 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 0.0 0.39269908169872414 0.7853981633974483 1.1780972450961724 1.5707963267948966 1.9634954084936207 2.356194490192345 2.748893571891069 3.141592653589793 0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586 + 0.0 0.0 0.0 1 From 94547f7c4b4128306a21a15d9adce1ac34d01c43 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 15:24:13 +0000 Subject: [PATCH 18/58] new type for centre_ --- include/openmc/mesh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index ea6b01aa60..8fdf62ecb6 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,7 +359,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - vector centre_; + array centre_; int set_grid(); From 15f45bcd1d7f6edf48c37213f0305803b41b584c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 15:29:10 +0000 Subject: [PATCH 19/58] substract before mapping --- src/mesh.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index c049549eb2..36e4d50490 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -976,6 +976,11 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( mapped_r[0] = std::hypot(r.x, r.y); mapped_r[2] = r[2]; + + mapped_r[0] += centre_[0]; + mapped_r[1] += centre_[1]; + mapped_r[2] += centre_[2]; + if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; } else { @@ -984,11 +989,6 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( mapped_r[1] += 2 * M_PI; } - // TODO: pass centre as argument - mapped_r[0] += centre_[0]; - mapped_r[1] += centre_[1]; - mapped_r[2] += centre_[2]; - MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_phi(idx[1]); @@ -1213,6 +1213,11 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( Position mapped_r; mapped_r[0] = r.norm(); + + mapped_r[0] += centre_[0]; + mapped_r[1] += centre_[1]; + mapped_r[2] += centre_[2]; + if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; mapped_r[2] = 0.0; @@ -1223,11 +1228,6 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[2] += 2 * M_PI; } - // TODO: pass centre as argument - mapped_r[0] += centre_[0]; - mapped_r[1] += centre_[1]; - mapped_r[2] += centre_[2]; - MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_theta(idx[1]); From 27e3bb427a3cf991dfa7f14680e8098835251443 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 15:36:58 +0000 Subject: [PATCH 20/58] array for origin CylindricalMesh --- include/openmc/mesh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 8fdf62ecb6..e6f5ce98b9 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -415,7 +415,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - vector centre_; + array centre_; int set_grid(); From 76fbdf7fbf627a177017853583c378192912d6b4 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 15:41:38 +0000 Subject: [PATCH 21/58] centre to origin --- include/openmc/mesh.h | 4 ++-- openmc/mesh.py | 52 +++++++++++++++++++++---------------------- src/mesh.cpp | 16 ++++++------- 3 files changed, 36 insertions(+), 36 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index e6f5ce98b9..9a0b5d5b84 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,7 +359,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array centre_; + array origin; int set_grid(); @@ -415,7 +415,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array centre_; + array origin_; int set_grid(); diff --git a/openmc/mesh.py b/openmc/mesh.py index f381d99b05..da00750c10 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,7 +1052,7 @@ class CylindricalMesh(StructuredMesh): self._r_grid = None self._phi_grid = [0.0, 2*pi] self._z_grid = None - self._centre = [0., 0., 0.] + self._origin = [0., 0., 0.] @property def dimension(self): @@ -1065,8 +1065,8 @@ class CylindricalMesh(StructuredMesh): return 3 @property - def centre(self): - return self._centre + def origin(self): + return self._origin @property def r_grid(self): @@ -1094,10 +1094,10 @@ class CylindricalMesh(StructuredMesh): for p in range(1, np + 1) for r in range(1, nr + 1)) - @centre.setter - def centre(self, coords): + @origin.setter + def origin(self, coords): cv.check_type('mesh r_grid', coords, Iterable, Real) - self._centre = np.asarray(coords) + self._origin = np.asarray(coords) @r_grid.setter def r_grid(self, grid): @@ -1240,8 +1240,8 @@ class CylindricalMesh(StructuredMesh): subelement = ET.SubElement(element, "z_grid") subelement.text = ' '.join(map(str, self.z_grid)) - subelement = ET.SubElement(element, "centre") - subelement.text = ' '.join(map(str, self.centre)) + subelement = ET.SubElement(element, "origin") + subelement.text = ' '.join(map(str, self.origin)) return element @@ -1266,7 +1266,7 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] - mesh.centre = [float(x) for x in get_text(elem, "centre").split()] + mesh.origin = [float(x) for x in get_text(elem, "origin").split()] return mesh @@ -1321,10 +1321,10 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) - # offset with centre - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2] + # offset with origin + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, @@ -1374,7 +1374,7 @@ class SphericalMesh(StructuredMesh): self._r_grid = None self._theta_grid = [0, pi] self._phi_grid = [0, 2*pi] - self._centre = [0., 0., 0.] + self._origin = [0., 0., 0.] @property def dimension(self): @@ -1387,8 +1387,8 @@ class SphericalMesh(StructuredMesh): return 3 @property - def centre(self): - return self._centre + def origin(self): + return self._origin @property def r_grid(self): @@ -1416,10 +1416,10 @@ class SphericalMesh(StructuredMesh): for t in range(1, nt + 1) for r in range(1, nr + 1)) - @centre.setter - def centre(self, coords): + @origin.setter + def origin(self, coords): cv.check_type('mesh r_grid', coords, Iterable, Real) - self._centre = np.asarray(coords) + self._origin = np.asarray(coords) @r_grid.setter def r_grid(self, grid): @@ -1492,8 +1492,8 @@ class SphericalMesh(StructuredMesh): subelement = ET.SubElement(element, "phi_grid") subelement.text = ' '.join(map(str, self.phi_grid)) - subelement = ET.SubElement(element, "centre") - subelement.text = ' '.join(map(str, self.centre)) + subelement = ET.SubElement(element, "origin") + subelement.text = ' '.join(map(str, self.origin)) return element @@ -1518,7 +1518,7 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] - mesh.centre = [float(x) for x in get_text(elem, "centre").split()] + mesh.origin = [float(x) for x in get_text(elem, "origin").split()] return mesh @@ -1575,10 +1575,10 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) pts_cartesian[:, 2] = r * np.cos(phi) - # offset with centre - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.centre[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.centre[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.centre[2] + # offset with origin + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, diff --git a/src/mesh.cpp b/src/mesh.cpp index 36e4d50490..982b2e2dc2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - centre_ = get_node_array(node, "centre"); + origin_ = get_node_array(node, "origin"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -977,9 +977,9 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( mapped_r[0] = std::hypot(r.x, r.y); mapped_r[2] = r[2]; - mapped_r[0] += centre_[0]; - mapped_r[1] += centre_[1]; - mapped_r[2] += centre_[2]; + mapped_r[0] += origin_[0]; + mapped_r[1] += origin_[1]; + mapped_r[2] += origin_[2]; if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - centre_ = get_node_array(node, "centre"); + origin_ = get_node_array(node, "origin"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1214,9 +1214,9 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[0] = r.norm(); - mapped_r[0] += centre_[0]; - mapped_r[1] += centre_[1]; - mapped_r[2] += centre_[2]; + mapped_r[0] += origin_[0]; + mapped_r[1] += origin_[1]; + mapped_r[2] += origin_[2]; if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; From b2d1729287c7174503c94e2b1cb5a39c2cb5a05e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 16:08:03 +0000 Subject: [PATCH 22/58] origin_ instead of origin --- include/openmc/mesh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 9a0b5d5b84..52f64701b5 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,7 +359,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array origin; + array origin_; int set_grid(); From adad6e7cd43a6d0a3e6ffd31f98bd7030158ee15 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 11 Oct 2022 11:21:31 -0500 Subject: [PATCH 23/58] Making origin type Positiong. Adding function to parse position from node. --- include/openmc/mesh.h | 5 ++--- include/openmc/xml_interface.h | 4 ++++ src/mesh.cpp | 4 ++-- src/xml_interface.cpp | 7 +++++++ 4 files changed, 15 insertions(+), 5 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 9a0b5d5b84..c690dc4fee 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -359,8 +359,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array origin; - + Position origin_; int set_grid(); @@ -415,7 +414,7 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - array origin_; + Position origin_; int set_grid(); diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h index d96bb6cc20..bd6554c134 100644 --- a/include/openmc/xml_interface.h +++ b/include/openmc/xml_interface.h @@ -9,6 +9,7 @@ #include "xtensor/xadapt.hpp" #include "xtensor/xarray.hpp" +#include "openmc/position.h" #include "openmc/vector.h" namespace openmc { @@ -49,5 +50,8 @@ xt::xarray get_node_xarray( return xt::adapt(v, shape); } +Position get_node_position( + pugi::xml_node node, const char* name, bool lowercase = false); + } // namespace openmc #endif // OPENMC_XML_INTERFACE_H diff --git a/src/mesh.cpp b/src/mesh.cpp index 982b2e2dc2..e9c185c8f1 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -955,7 +955,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); - origin_ = get_node_array(node, "origin"); + origin_ = get_node_position(node, "origin"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1193,7 +1193,7 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); - origin_ = get_node_array(node, "origin"); + origin_ = get_node_position(node, "origin"); if (int err = set_grid()) { fatal_error(openmc_err_msg); diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index cf3b981e51..6ce465bafb 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -48,4 +48,11 @@ bool get_node_value_bool(pugi::xml_node node, const char* name) return false; } +Position get_node_position( + pugi::xml_node node, const char* name, bool lowercase) +{ + vector arr = get_node_array(node, name, lowercase); + return Position(arr); +} + } // namespace openmc From a97b202c2d5113795ff39ab95827fd61ecf3ba7f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 16:41:43 +0000 Subject: [PATCH 24/58] updated to_hdf5 --- src/mesh.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index e9c185c8f1..57fe49ba6b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1178,6 +1178,7 @@ void CylindricalMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "r_grid", grid_[0]); write_dataset(mesh_group, "phi_grid", grid_[1]); write_dataset(mesh_group, "z_grid", grid_[2]); + write_dataset(mesh_group, "origin", origin_); close_group(mesh_group); } @@ -1447,6 +1448,7 @@ void SphericalMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "r_grid", grid_[0]); write_dataset(mesh_group, "theta_grid", grid_[1]); write_dataset(mesh_group, "phi_grid", grid_[2]); + write_dataset(mesh_group, "origin", origin_); close_group(mesh_group); } From 3792c6fb70786638d39c88201fc248464c0bd905 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 17:11:54 +0000 Subject: [PATCH 25/58] fixed unit test --- tests/unit_tests/test_spherical_mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index cdeb188ece..4d389a6bde 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -1,8 +1,8 @@ import openmc import numpy as np -def test_centre_read_write_to_xml(): - """Tests that the centre attribute can be written and read back to XML +def test_origin_read_write_to_xml(): + """Tests that the origin attribute can be written and read back to XML """ # build mesh = openmc.SphericalMesh() @@ -10,7 +10,7 @@ def test_centre_read_write_to_xml(): mesh.theta_grid = [1, 2, 3] mesh.r_grid = [1, 2, 3] - mesh.centre = [0.1, 0.2, 0.3] + mesh.origin = [0.1, 0.2, 0.3] tally = openmc.Tally() @@ -27,4 +27,4 @@ def test_centre_read_write_to_xml(): new_tallies = openmc.Tallies.from_xml() new_tally = new_tallies[0] new_mesh = new_tally.filters[0].mesh - assert np.allclose(new_mesh.centre, mesh.centre) \ No newline at end of file + assert np.allclose(new_mesh.origin, mesh.origin) \ No newline at end of file From 48d75accdd3d8ddd80e560c56b64a49e3c7b21dd Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 17:12:55 +0000 Subject: [PATCH 26/58] fixed cylindrical mesh test --- tests/unit_tests/test_cylindrical_mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 89322d7045..952e8cf5a2 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -1,8 +1,8 @@ import openmc import numpy as np -def test_centre_read_write_to_xml(): - """Tests that the centre attribute can be written and read back to XML +def test_origin_read_write_to_xml(): + """Tests that the origin attribute can be written and read back to XML """ # build mesh = openmc.CylindricalMesh() @@ -10,7 +10,7 @@ def test_centre_read_write_to_xml(): mesh.z_grid = [1, 2, 3] mesh.r_grid = [1, 2, 3] - mesh.centre = [0.1, 0.2, 0.3] + mesh.origin = [0.1, 0.2, 0.3] tally = openmc.Tally() @@ -27,4 +27,4 @@ def test_centre_read_write_to_xml(): new_tallies = openmc.Tallies.from_xml() new_tally = new_tallies[0] new_mesh = new_tally.filters[0].mesh - assert np.allclose(new_mesh.centre, mesh.centre) \ No newline at end of file + assert np.allclose(new_mesh.origin, mesh.origin) \ No newline at end of file From c6b89feaef6fc8093fc012e0a4a426541c7d434d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 11 Oct 2022 17:14:21 +0000 Subject: [PATCH 27/58] fixed gold standard --- tests/regression_tests/filter_mesh/inputs_true.dat | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index c42f3cd439..8b08a2a3c1 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -57,13 +57,13 @@ 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586 -7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5 - 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 0.0 0.39269908169872414 0.7853981633974483 1.1780972450961724 1.5707963267948966 1.9634954084936207 2.356194490192345 2.748893571891069 3.141592653589793 0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586 - 0.0 0.0 0.0 + 0.0 0.0 0.0 1 From 444a80012bf6d246842f133a89943ccf1004fd87 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 11 Oct 2022 23:38:13 -0500 Subject: [PATCH 28/58] Naming tally in activation unit test --- tests/unit_tests/test_deplete_activation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 1842ad8ac8..2dc911c86f 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -33,7 +33,7 @@ def model(): ) model.settings.run_mode = 'fixed source' - rx_tally = openmc.Tally() + rx_tally = openmc.Tally(name='activation tally') rx_tally.scores = ['(n,gamma)'] model.tallies.append(rx_tally) @@ -53,7 +53,8 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts # Determine (n.gamma) reaction rate using initial run sp = model.run() with openmc.StatePoint(sp) as sp: - tally = sp.tallies[1] + print(sp.tallies) + tally = sp.get_tally(name='activation tally') capture_rate = tally.mean.flat[0] # Create one-nuclide depletion chain From b0743dfdea203372844e3d4901d8d0f8b2150b0c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 21 Oct 2022 11:56:21 +0200 Subject: [PATCH 29/58] fixed cv.check_type origin --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index da00750c10..8078869889 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1096,7 +1096,7 @@ class CylindricalMesh(StructuredMesh): @origin.setter def origin(self, coords): - cv.check_type('mesh r_grid', coords, Iterable, Real) + cv.check_type('mesh origin', coords, Iterable, Real) self._origin = np.asarray(coords) @r_grid.setter @@ -1418,7 +1418,7 @@ class SphericalMesh(StructuredMesh): @origin.setter def origin(self, coords): - cv.check_type('mesh r_grid', coords, Iterable, Real) + cv.check_type('mesh origin', coords, Iterable, Real) self._origin = np.asarray(coords) @r_grid.setter From b874d83634962c9e1fcae2d0ffc360437e11abf5 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 21 Oct 2022 15:55:08 +0200 Subject: [PATCH 30/58] only translate when origin is not 0, 0, 0 --- openmc/mesh.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 8078869889..b4cd04d511 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1322,9 +1322,10 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) # offset with origin - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] + if self.origin != [0.0, 0.0, 0.0]: + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, @@ -1576,9 +1577,10 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 2] = r * np.cos(phi) # offset with origin - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] + if self.origin != [0.0, 0.0, 0.0]: + pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] + pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] + pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, From f60f95feb4e49d4b46c6aa5ec63b7ed4ddeeec87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Fri, 18 Nov 2022 09:03:35 -0500 Subject: [PATCH 31/58] [skip ci] removed # TODOs --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index b4cd04d511..79ab7fd5a6 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -165,7 +165,7 @@ class StructuredMesh(MeshBase): pass @property - def vertices(self): # TODO needs changing + def vertices(self): """Return coordinates of mesh vertices. Returns @@ -178,7 +178,7 @@ class StructuredMesh(MeshBase): return np.stack(np.meshgrid(*self._grids, indexing='ij'), axis=-1) @property - def centroids(self): # TODO needs changing + def centroids(self): """Return coordinates of mesh element centroids. Returns From c5affa805e06ce46276dc1c148888b3540b1e29d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 21 Nov 2022 09:37:57 -0500 Subject: [PATCH 32/58] Review suggestions Co-authored-by: Patrick Shriwise --- src/mesh.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 57fe49ba6b..3849f4d031 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -977,9 +977,7 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( mapped_r[0] = std::hypot(r.x, r.y); mapped_r[2] = r[2]; - mapped_r[0] += origin_[0]; - mapped_r[1] += origin_[1]; - mapped_r[2] += origin_[2]; + mapped_r -= origin_; if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; From 6d787d5fb0e2ed7790e83f2ea8d00560f60c292e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 21 Nov 2022 09:39:01 -0500 Subject: [PATCH 33/58] Apply suggestions from code review Co-authored-by: Patrick Shriwise --- openmc/mesh.py | 6 +++--- src/mesh.cpp | 4 +--- tests/unit_tests/test_deplete_activation.py | 1 - 3 files changed, 4 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 79ab7fd5a6..3936572f9d 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1322,7 +1322,7 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) # offset with origin - if self.origin != [0.0, 0.0, 0.0]: + if any(self.origin != (0.0, 0.0, 0.0)): pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] @@ -1375,7 +1375,7 @@ class SphericalMesh(StructuredMesh): self._r_grid = None self._theta_grid = [0, pi] self._phi_grid = [0, 2*pi] - self._origin = [0., 0., 0.] + self._origin = (0., 0., 0.) @property def dimension(self): @@ -1577,7 +1577,7 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 2] = r * np.cos(phi) # offset with origin - if self.origin != [0.0, 0.0, 0.0]: + if any(self.origin != (0.0, 0.0, 0.0)): pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] diff --git a/src/mesh.cpp b/src/mesh.cpp index 3849f4d031..82644db08b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1213,9 +1213,7 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( mapped_r[0] = r.norm(); - mapped_r[0] += origin_[0]; - mapped_r[1] += origin_[1]; - mapped_r[2] += origin_[2]; + mapped_r -= origin_; if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 2dc911c86f..64167c8af9 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -53,7 +53,6 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts # Determine (n.gamma) reaction rate using initial run sp = model.run() with openmc.StatePoint(sp) as sp: - print(sp.tallies) tally = sp.get_tally(name='activation tally') capture_rate = tally.mean.flat[0] From 8681751238a27d4b94a75c446903c035de6f72dd Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 09:45:11 -0500 Subject: [PATCH 34/58] check length origin --- openmc/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3936572f9d..0eba6871b1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1097,6 +1097,7 @@ class CylindricalMesh(StructuredMesh): @origin.setter def origin(self, coords): cv.check_type('mesh origin', coords, Iterable, Real) + cv.check_length("mesh origin", coords, 3) self._origin = np.asarray(coords) @r_grid.setter @@ -1420,6 +1421,7 @@ class SphericalMesh(StructuredMesh): @origin.setter def origin(self, coords): cv.check_type('mesh origin', coords, Iterable, Real) + cv.check_length("mesh origin", coords, 3) self._origin = np.asarray(coords) @r_grid.setter From df1a45b65407f0dc67e286600aac58d1b8744682 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 09:46:13 -0500 Subject: [PATCH 35/58] origin from_hdf5 --- openmc/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 0eba6871b1..64217fa777 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1145,6 +1145,7 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = group['r_grid'][()] mesh.phi_grid = group['phi_grid'][()] mesh.z_grid = group['z_grid'][()] + mesh.origin = group['origin'][()] return mesh @@ -1469,6 +1470,7 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = group['r_grid'][()] mesh.theta_grid = group['theta_grid'][()] mesh.phi_grid = group['phi_grid'][()] + mesh.origin = group['origin'][()] return mesh From b63ec9cdc3ee988afe1051126cdba2b7d83ff629 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 09:48:42 -0500 Subject: [PATCH 36/58] read mesh optional --- openmc/mesh.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 64217fa777..0d92208aa1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1145,7 +1145,8 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = group['r_grid'][()] mesh.phi_grid = group['phi_grid'][()] mesh.z_grid = group['z_grid'][()] - mesh.origin = group['origin'][()] + if 'origin' in group: + mesh.origin = group['origin'][()] return mesh @@ -1268,7 +1269,8 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] - mesh.origin = [float(x) for x in get_text(elem, "origin").split()] + if "origin" in elem.attrib: + mesh.origin = [float(x) for x in get_text(elem, "origin").split()] return mesh @@ -1470,7 +1472,8 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = group['r_grid'][()] mesh.theta_grid = group['theta_grid'][()] mesh.phi_grid = group['phi_grid'][()] - mesh.origin = group['origin'][()] + if 'origin' in group: + mesh.origin = group['origin'][()] return mesh @@ -1523,7 +1526,8 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] - mesh.origin = [float(x) for x in get_text(elem, "origin").split()] + if "origin" in elem.attrib: + mesh.origin = [float(x) for x in get_text(elem, "origin").split()] return mesh From f51f18176135926ea1c4be5f2287c95a8741c8fd Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 11:53:49 -0500 Subject: [PATCH 37/58] fixed optional xml reading --- openmc/mesh.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 0d92208aa1..31b7789ee1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1269,8 +1269,7 @@ class CylindricalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()] - if "origin" in elem.attrib: - mesh.origin = [float(x) for x in get_text(elem, "origin").split()] + mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()] return mesh @@ -1526,8 +1525,7 @@ class SphericalMesh(StructuredMesh): mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()] mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()] mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()] - if "origin" in elem.attrib: - mesh.origin = [float(x) for x in get_text(elem, "origin").split()] + mesh.origin = [float(x) for x in get_text(elem, "origin", default=[0., 0., 0.]).split()] return mesh From 3325667c51538748e2078261fb9457c769b11e3e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 12:22:39 -0500 Subject: [PATCH 38/58] minor version bump for statepoint files --- include/openmc/constants.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index fa2251b84d..a97180fafa 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -24,7 +24,7 @@ using double_4dvec = vector>>>; constexpr int HDF5_VERSION[] {3, 0}; // Version numbers for binary files -constexpr array VERSION_STATEPOINT {18, 0}; +constexpr array VERSION_STATEPOINT {18, 1}; constexpr array VERSION_PARTICLE_RESTART {2, 0}; constexpr array VERSION_TRACK {3, 0}; constexpr array VERSION_SUMMARY {6, 0}; From 10a95186e391517ff0f5d1f7487d6aa4c4cd8d93 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 21 Nov 2022 14:07:45 -0500 Subject: [PATCH 39/58] origin is tuple + fixed any --- openmc/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 31b7789ee1..615514c683 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,7 +1052,7 @@ class CylindricalMesh(StructuredMesh): self._r_grid = None self._phi_grid = [0.0, 2*pi] self._z_grid = None - self._origin = [0., 0., 0.] + self._origin = (0., 0., 0.) @property def dimension(self): @@ -1325,7 +1325,7 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) # offset with origin - if any(self.origin != (0.0, 0.0, 0.0)): + if any([coord != 0 for coord in self.origin]): pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] @@ -1583,7 +1583,7 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 2] = r * np.cos(phi) # offset with origin - if any(self.origin != (0.0, 0.0, 0.0)): + if any([coord != 0 for coord in self.origin]): pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] From f0721ac5351eec9425ecf4e33a90b5ac13eba1cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 21 Nov 2022 14:30:36 -0500 Subject: [PATCH 40/58] Apply suggestions from code review Co-authored-by: Patrick Shriwise --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 615514c683..307d9d952c 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1052,7 +1052,7 @@ class CylindricalMesh(StructuredMesh): self._r_grid = None self._phi_grid = [0.0, 2*pi] self._z_grid = None - self._origin = (0., 0., 0.) + self.origin = (0., 0., 0.) @property def dimension(self): @@ -1378,7 +1378,7 @@ class SphericalMesh(StructuredMesh): self._r_grid = None self._theta_grid = [0, pi] self._phi_grid = [0, 2*pi] - self._origin = (0., 0., 0.) + self.origin = (0., 0., 0.) @property def dimension(self): From 11eee756dd953d2490fa8b2975341e7f64328f79 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 26 Jan 2023 14:13:39 -0600 Subject: [PATCH 41/58] Corrections to application of origin for sphere and cyl mesh --- include/openmc/mesh.h | 3 +-- src/mesh.cpp | 27 +++++++++++++++------------ 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index c690dc4fee..9554a51436 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -241,6 +241,7 @@ public: xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh std::array shape_; //!< Number of mesh elements in each dimension + Position origin_ {0.0, 0.0, 0.0}; protected: }; @@ -359,7 +360,6 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - Position origin_; int set_grid(); @@ -414,7 +414,6 @@ public: void to_hdf5(hid_t group) const override; array, 3> grid_; - Position origin_; int set_grid(); diff --git a/src/mesh.cpp b/src/mesh.cpp index 82644db08b..e0fa305edd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -412,7 +412,6 @@ template void StructuredMesh::raytrace_mesh( Position r0, Position r1, const Direction& u, T tally) const { - // TODO: when c++-17 is available, use "if constexpr ()" to compile-time // enable/disable tally calls for now, T template type needs to provide both // surface and track methods, which might be empty. modern optimizing @@ -445,6 +444,11 @@ void StructuredMesh::raytrace_mesh( return; } + // translate start and end positions, + // this needs to come after the get_indices call because it does its own translation + r0 -= origin_; + r1 -= origin_; + // Calculate initial distances to next surfaces in all three dimensions std::array distances; for (int k = 0; k < n; ++k) { @@ -972,13 +976,12 @@ std::string CylindricalMesh::get_mesh_type() const StructuredMesh::MeshIndex CylindricalMesh::get_indices( Position r, bool& in_mesh) const { - Position mapped_r; + r -= origin_; + Position mapped_r; mapped_r[0] = std::hypot(r.x, r.y); mapped_r[2] = r[2]; - mapped_r -= origin_; - if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; } else { @@ -1090,22 +1093,23 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary( const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const { + Position r = r0 - origin_; if (i == 0) { return std::min( - MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])), - MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1))); + MeshDistance(ijk[i] + 1, true, find_r_crossing(r, u, l, ijk[i])), + MeshDistance(ijk[i] - 1, false, find_r_crossing(r, u, l, ijk[i] - 1))); } else if (i == 1) { return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true, - find_phi_crossing(r0, u, l, ijk[i])), + find_phi_crossing(r, u, l, ijk[i])), MeshDistance(sanitize_phi(ijk[i] - 1), false, - find_phi_crossing(r0, u, l, ijk[i] - 1))); + find_phi_crossing(r, u, l, ijk[i] - 1))); } else { - return find_z_crossing(r0, u, l, ijk[i]); + return find_z_crossing(r, u, l, ijk[i]); } } @@ -1209,12 +1213,11 @@ std::string SphericalMesh::get_mesh_type() const StructuredMesh::MeshIndex SphericalMesh::get_indices( Position r, bool& in_mesh) const { + r -= origin_; + Position mapped_r; - mapped_r[0] = r.norm(); - mapped_r -= origin_; - if (mapped_r[0] < FP_PRECISION) { mapped_r[1] = 0.0; mapped_r[2] = 0.0; From b22c7e8353781552d1cedddd05b83ad2543c5714 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 26 Jan 2023 14:14:06 -0600 Subject: [PATCH 42/58] Adding test cases for cylindrical mesh --- tests/unit_tests/test_cylindrical_mesh.py | 88 ++++++++++++++++++++--- 1 file changed, 78 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 952e8cf5a2..e67f21b15c 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -1,16 +1,42 @@ +from itertools import product, permutations + import openmc import numpy as np -def test_origin_read_write_to_xml(): - """Tests that the origin attribute can be written and read back to XML - """ +import pytest + +geom_size = 5 + +@pytest.fixture() +def model(): + openmc.reset_auto_ids() + + water = openmc.Material(name='water') + water.add_element('H', 2.0) + water.add_element('O', 1.0) + water.set_density('g/cc', 1.0) + + rpp = openmc.model.RectangularParallelepiped(*([-geom_size, geom_size] * 3), + boundary_type='vacuum') + + cell = openmc.Cell(region=-rpp, fill=water) + + geom = openmc.Geometry([cell]) + + source = openmc.Source() + source.space = openmc.stats.Point() + source.energy = openmc.stats.Discrete([10000], [1.0]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 10 + settings.run_mode = 'fixed source' + # build mesh = openmc.CylindricalMesh() - mesh.phi_grid = [1, 2, 3] - mesh.z_grid = [1, 2, 3] - mesh.r_grid = [1, 2, 3] - - mesh.origin = [0.1, 0.2, 0.3] + mesh.phi_grid = np.linspace(0, 2*np.pi, 21) + mesh.z_grid = np.linspace(-geom_size, geom_size, 11) + mesh.r_grid = np.linspace(0, geom_size, geom_size) tally = openmc.Tally() @@ -21,10 +47,52 @@ def test_origin_read_write_to_xml(): tallies = openmc.Tallies([tally]) - tallies.export_to_xml() + return openmc.Model(geometry=geom, settings=settings, tallies=tallies) + +def test_origin_read_write_to_xml(run_in_tmpdir, model): + """Tests that the origin attribute can be written and read back to XML + """ + mesh = model.tallies[0].filters[0].mesh + mesh.origin = [0.1, 0.2, 0.3] + model.tallies.export_to_xml() # read back new_tallies = openmc.Tallies.from_xml() new_tally = new_tallies[0] new_mesh = new_tally.filters[0].mesh - assert np.allclose(new_mesh.origin, mesh.origin) \ No newline at end of file + np.testing.assert_equal(new_mesh.origin, mesh.origin) + +estimators = ('tracklength', 'collision') +origins = permutations([-geom_size, 0, 0]) + +test_cases = product(estimators, origins) + +def label(p): + if isinstance(p, tuple): + return f'origin:{p}' + if isinstance(p, str): + return f'estimator:{p}' + +@pytest.mark.parametrize('estimator,origin', test_cases, ids=label) +def test_offset_mesh(model, estimator, origin): + """Tests that the mesh has been moved based on tally results + """ + mesh = model.tallies[0].filters[0].mesh + model.tallies[0].estimator = estimator + # move the center of the cylinder mesh upwards + mesh.origin = origin + + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + tally = sp.tallies[1] + + # we've translated half of the cylinder mesh above the model, + # so ensure that half of the bins are populated + assert np.count_nonzero(tally.mean) == tally.mean.size / 2 + + # check that the lower half of the mesh contains a tally result + # and that the upper half is zero + mean = tally.mean.reshape(mesh.dimension[::-1]).T + # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 + # assert np.count_nonzero(mean[:, :, 5:]) == 0 From 3703fdbdd746675203e1c5477a29b521cdb9a827 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 26 Jan 2023 14:15:08 -0600 Subject: [PATCH 43/58] Adding origin to sphere and cyl mesh str --- openmc/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 307d9d952c..47d34ebe2e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1119,6 +1119,7 @@ class CylindricalMesh(StructuredMesh): fmt = '{0: <16}{1}{2}\n' string = super().__repr__() string += fmt.format('\tDimensions', '=\t', self.n_dimension) + string += fmt.format('\tOrigin', '=\t', self.origin) r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid) string += fmt.format('\tN R pnts:', '=\t', r_grid_str) if self._r_grid is not None: @@ -1445,6 +1446,7 @@ class SphericalMesh(StructuredMesh): fmt = '{0: <16}{1}{2}\n' string = super().__repr__() string += fmt.format('\tDimensions', '=\t', self.n_dimension) + string += fmt.format('\tOrigin', '=\t', self.origin) r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid) string += fmt.format('\tN R pnts:', '=\t', r_grid_str) if self._r_grid is not None: From 17d88d8c845bb84cfbda4f674ea0919d64b79154 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 30 Jan 2023 13:57:44 -0500 Subject: [PATCH 44/58] added tests for spherical mesh --- tests/unit_tests/test_spherical_mesh.py | 88 ++++++++++++++++++++++--- 1 file changed, 78 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 4d389a6bde..711570fb31 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -1,16 +1,42 @@ +from itertools import product, permutations + import openmc import numpy as np -def test_origin_read_write_to_xml(): - """Tests that the origin attribute can be written and read back to XML - """ +import pytest + +geom_size = 5 + +@pytest.fixture() +def model(): + openmc.reset_auto_ids() + + water = openmc.Material(name='water') + water.add_element('H', 2.0) + water.add_element('O', 1.0) + water.set_density('g/cc', 1.0) + + rpp = openmc.model.RectangularParallelepiped(*([-geom_size, geom_size] * 3), + boundary_type='vacuum') + + cell = openmc.Cell(region=-rpp, fill=water) + + geom = openmc.Geometry([cell]) + + source = openmc.Source() + source.space = openmc.stats.Point() + source.energy = openmc.stats.Discrete([10000], [1.0]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 10 + settings.run_mode = 'fixed source' + # build mesh = openmc.SphericalMesh() - mesh.phi_grid = [1, 2, 3] - mesh.theta_grid = [1, 2, 3] - mesh.r_grid = [1, 2, 3] - - mesh.origin = [0.1, 0.2, 0.3] + mesh.phi_grid = np.linspace(0, 2*np.pi, 21) + mesh.theta_grid = np.linspace(0, 2*np.pi, 11) + mesh.r_grid = np.linspace(0, geom_size, geom_size) tally = openmc.Tally() @@ -21,10 +47,52 @@ def test_origin_read_write_to_xml(): tallies = openmc.Tallies([tally]) - tallies.export_to_xml() + return openmc.Model(geometry=geom, settings=settings, tallies=tallies) + +def test_origin_read_write_to_xml(run_in_tmpdir, model): + """Tests that the origin attribute can be written and read back to XML + """ + mesh = model.tallies[0].filters[0].mesh + mesh.origin = [0.1, 0.2, 0.3] + model.tallies.export_to_xml() # read back new_tallies = openmc.Tallies.from_xml() new_tally = new_tallies[0] new_mesh = new_tally.filters[0].mesh - assert np.allclose(new_mesh.origin, mesh.origin) \ No newline at end of file + np.testing.assert_equal(new_mesh.origin, mesh.origin) + +estimators = ('tracklength', 'collision') +origins = permutations([-geom_size, 0, 0]) + +test_cases = product(estimators, origins) + +def label(p): + if isinstance(p, tuple): + return f'origin:{p}' + if isinstance(p, str): + return f'estimator:{p}' + +@pytest.mark.parametrize('estimator,origin', test_cases, ids=label) +def test_offset_mesh(model, estimator, origin): + """Tests that the mesh has been moved based on tally results + """ + mesh = model.tallies[0].filters[0].mesh + model.tallies[0].estimator = estimator + # move the center of the spherical mesh upwards + mesh.origin = origin + + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + tally = sp.tallies[1] + + # we've translated half of the spherical mesh above the model, + # so ensure that half of the bins are populated + assert np.count_nonzero(tally.mean) == tally.mean.size / 2 + + # check that the lower half of the mesh contains a tally result + # and that the upper half is zero + mean = tally.mean.reshape(mesh.dimension[::-1]).T + # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 + # assert np.count_nonzero(mean[:, :, 5:]) == 0 From 1e94a0c7d7126dc68e0dfdda933a9142e8cd1f42 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 2 Feb 2023 09:39:23 -0500 Subject: [PATCH 45/58] fixed theta_grid spherical mesh --- tests/unit_tests/test_spherical_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 711570fb31..d6ac435c2e 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -35,7 +35,7 @@ def model(): # build mesh = openmc.SphericalMesh() mesh.phi_grid = np.linspace(0, 2*np.pi, 21) - mesh.theta_grid = np.linspace(0, 2*np.pi, 11) + mesh.theta_grid = np.linspace(0, np.pi, 11) mesh.r_grid = np.linspace(0, geom_size, geom_size) tally = openmc.Tally() From af16239684024712beec6ace9b74e3024ad34b20 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 8 Feb 2023 10:17:59 -0600 Subject: [PATCH 46/58] Updating mesh centre tests. --- tests/unit_tests/test_cylindrical_mesh.py | 21 ++++++++++----- tests/unit_tests/test_spherical_mesh.py | 31 +++++++++++++++-------- 2 files changed, 35 insertions(+), 17 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index e67f21b15c..a9d3fd56a8 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -28,7 +28,7 @@ def model(): source.energy = openmc.stats.Discrete([10000], [1.0]) settings = openmc.Settings() - settings.particles = 1000 + settings.particles = 2000 settings.batches = 10 settings.run_mode = 'fixed source' @@ -43,7 +43,7 @@ def model(): mesh_filter = openmc.MeshFilter(mesh) tally.filters.append(mesh_filter) - tally.scores.append("heating") + tally.scores.append("flux") tallies = openmc.Tallies([tally]) @@ -63,7 +63,8 @@ def test_origin_read_write_to_xml(run_in_tmpdir, model): np.testing.assert_equal(new_mesh.origin, mesh.origin) estimators = ('tracklength', 'collision') -origins = permutations([-geom_size, 0, 0]) +origins = set(permutations((-geom_size, 0, 0))) +origins |= set(permutations((geom_size, 0, 0))) test_cases = product(estimators, origins) @@ -91,8 +92,14 @@ def test_offset_mesh(model, estimator, origin): # so ensure that half of the bins are populated assert np.count_nonzero(tally.mean) == tally.mean.size / 2 - # check that the lower half of the mesh contains a tally result - # and that the upper half is zero + # check that the half of the mesh that is outside of the geometry + # contains the zero values mean = tally.mean.reshape(mesh.dimension[::-1]).T - # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 - # assert np.count_nonzero(mean[:, :, 5:]) == 0 + centroids = mesh.centroids + for ijk in mesh.indices: + i, j, k = np.array(ijk) - 1 + print(centroids[:, i, j, k]) + if model.geometry.find(centroids[:, i, j, k]): + mean[i, j, k] == 0.0 + else: + mean[i, j, k] != 0.0 diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index d6ac435c2e..3cd8bbc5d7 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -28,14 +28,14 @@ def model(): source.energy = openmc.stats.Discrete([10000], [1.0]) settings = openmc.Settings() - settings.particles = 1000 + settings.particles = 2000 settings.batches = 10 settings.run_mode = 'fixed source' # build mesh = openmc.SphericalMesh() - mesh.phi_grid = np.linspace(0, 2*np.pi, 21) - mesh.theta_grid = np.linspace(0, np.pi, 11) + mesh.phi_grid = np.linspace(0, 2*np.pi, 13) + mesh.theta_grid = np.linspace(0, np.pi, 7) mesh.r_grid = np.linspace(0, geom_size, geom_size) tally = openmc.Tally() @@ -43,7 +43,7 @@ def model(): mesh_filter = openmc.MeshFilter(mesh) tally.filters.append(mesh_filter) - tally.scores.append("heating") + tally.scores.append("flux") tallies = openmc.Tallies([tally]) @@ -63,7 +63,12 @@ def test_origin_read_write_to_xml(run_in_tmpdir, model): np.testing.assert_equal(new_mesh.origin, mesh.origin) estimators = ('tracklength', 'collision') -origins = permutations([-geom_size, 0, 0]) +# TODO: determine why this is needed for spherical mesh +# but not cylindrical mesh +offset = geom_size + 0.001 + +origins = set(permutations((-offset, 0, 0))) +origins |= set(permutations((offset, 0, 0))) test_cases = product(estimators, origins) @@ -79,7 +84,7 @@ def test_offset_mesh(model, estimator, origin): """ mesh = model.tallies[0].filters[0].mesh model.tallies[0].estimator = estimator - # move the center of the spherical mesh upwards + # move the center of the spherical mesh mesh.origin = origin sp_filename = model.run() @@ -91,8 +96,14 @@ def test_offset_mesh(model, estimator, origin): # so ensure that half of the bins are populated assert np.count_nonzero(tally.mean) == tally.mean.size / 2 - # check that the lower half of the mesh contains a tally result - # and that the upper half is zero + # check that the half of the mesh that is outside of the geometry + # contains the zero values mean = tally.mean.reshape(mesh.dimension[::-1]).T - # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 - # assert np.count_nonzero(mean[:, :, 5:]) == 0 + centroids = mesh.centroids + for ijk in mesh.indices: + i, j, k = np.array(ijk) - 1 + print(centroids[:, i, j, k]) + if model.geometry.find(centroids[:, i, j, k]): + mean[i, j, k] == 0.0 + else: + mean[i, j, k] != 0.0 From c634a1467445352974e391cdfe7e80d936d3c6a9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 17 Feb 2023 10:56:42 -0500 Subject: [PATCH 47/58] added origin to docstrings --- openmc/mesh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 10d07000e5..06f0246768 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1049,6 +1049,8 @@ class CylindricalMesh(StructuredMesh): The default value is [0, 2π], i.e. the full phi range. z_grid : numpy.ndarray 1-D array of mesh boundary points along the z-axis. + origin : tuple + The (x,y,z) origin of the mesh in cartesian coordinates indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] @@ -1376,6 +1378,8 @@ class SphericalMesh(StructuredMesh): phi_grid : numpy.ndarray 1-D array of mesh boundary points along the phi-axis in radians. The default value is [0, 2π], i.e. the full phi range. + origin : tuple + The (x,y,z) origin of the mesh in cartesian coordinates indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] From f961fc90c2bc5c0565be5179a6a4c789cc553729 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Tue, 21 Feb 2023 17:42:06 -0500 Subject: [PATCH 48/58] [skip ci] Modified the XML specs --- docs/source/io_formats/tallies.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 8594bb11c3..897941ad40 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -351,6 +351,9 @@ attributes/sub-elements: :theta_grid: The mesh divisions along the theta-axis. (For spherical mesh only.) + :origin: + The origin in cartesian coordinates. (For cylindrical and spherical meshes only.) + :library: The mesh library used to represent an unstructured mesh. This can be either "moab" or "libmesh". (For unstructured mesh only.) From 2f9f500542fef29ddf5ec1b408e9fd1e10c37f03 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Feb 2023 16:10:27 -0500 Subject: [PATCH 49/58] removed check and always offset origin --- openmc/mesh.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 06f0246768..2c74c7bfad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1337,10 +1337,9 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) # offset with origin - if any([coord != 0 for coord in self.origin]): - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] + pts_cartesian[:, 0] += pts_cartesian[:, 0] + pts_cartesian[:, 1] += pts_cartesian[:, 1] + pts_cartesian[:, 2] += pts_cartesian[:, 2] return super().write_data_to_vtk( points=pts_cartesian, @@ -1598,10 +1597,9 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 2] = r * np.cos(phi) # offset with origin - if any([coord != 0 for coord in self.origin]): - pts_cartesian[:, 0] = pts_cartesian[:, 0] + self.origin[0] - pts_cartesian[:, 1] = pts_cartesian[:, 1] + self.origin[1] - pts_cartesian[:, 2] = pts_cartesian[:, 2] + self.origin[2] + pts_cartesian[:, 0] += self.origin[0] + pts_cartesian[:, 1] += self.origin[1] + pts_cartesian[:, 2] += self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, From c61e735a706b3bd251010ec8fa888cda30ebb596 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Feb 2023 16:10:55 -0500 Subject: [PATCH 50/58] empty_like instead of copy --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2c74c7bfad..9bed9200ce 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1331,7 +1331,7 @@ class CylindricalMesh(StructuredMesh): for r in self.r_grid ] ) - pts_cartesian = np.copy(pts_cylindrical) + pts_cartesian = np.empty_like(pts_cylindrical) r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) @@ -1590,7 +1590,7 @@ class SphericalMesh(StructuredMesh): for r in self.r_grid ] ) - pts_cartesian = np.copy(pts_spherical) + pts_cartesian = np.empty_like(pts_spherical) r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) From 9ca6f1b65b469938b5bb067c1c6d9d11a21f3623 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Feb 2023 16:13:31 -0500 Subject: [PATCH 51/58] fixed docstrings for origin --- openmc/mesh.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 9bed9200ce..db57a3a9ec 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1049,8 +1049,9 @@ class CylindricalMesh(StructuredMesh): The default value is [0, 2π], i.e. the full phi range. z_grid : numpy.ndarray 1-D array of mesh boundary points along the z-axis. - origin : tuple - The (x,y,z) origin of the mesh in cartesian coordinates + origin : numpy.ndarray + 1-D array of length 3 the (x,y,z) origin of the mesh in + cartesian coordinates indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] @@ -1377,8 +1378,9 @@ class SphericalMesh(StructuredMesh): phi_grid : numpy.ndarray 1-D array of mesh boundary points along the phi-axis in radians. The default value is [0, 2π], i.e. the full phi range. - origin : tuple - The (x,y,z) origin of the mesh in cartesian coordinates + origin : numpy.ndarray + 1-D array of length 3 the (x,y,z) origin of the mesh in + cartesian coordinates indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] From 972e1c12ec4cc2001202409bf8227454f8c83abb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 27 Feb 2023 16:22:10 -0500 Subject: [PATCH 52/58] [skip ci] Minor refactoring Co-authored-by: Ethan Peterson --- openmc/mesh.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index db57a3a9ec..a7c2eb1729 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1594,14 +1594,9 @@ class SphericalMesh(StructuredMesh): ) pts_cartesian = np.empty_like(pts_spherical) r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] - pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) - pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) - pts_cartesian[:, 2] = r * np.cos(phi) - - # offset with origin - pts_cartesian[:, 0] += self.origin[0] - pts_cartesian[:, 1] += self.origin[1] - pts_cartesian[:, 2] += self.origin[2] + pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) + self.origin[0] + pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1] + pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2] return super().write_data_to_vtk( points=pts_cartesian, From 308be3f1eac95143eec2f3f40c36adb6dfa351c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Mon, 27 Feb 2023 16:23:45 -0500 Subject: [PATCH 53/58] Update openmc/mesh.py Co-authored-by: Ethan Peterson --- openmc/mesh.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index a7c2eb1729..039f302e92 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1334,14 +1334,10 @@ class CylindricalMesh(StructuredMesh): ) pts_cartesian = np.empty_like(pts_cylindrical) r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] - pts_cartesian[:, 0] = r * np.cos(phi) - pts_cartesian[:, 1] = r * np.sin(phi) - - # offset with origin - pts_cartesian[:, 0] += pts_cartesian[:, 0] - pts_cartesian[:, 1] += pts_cartesian[:, 1] - pts_cartesian[:, 2] += pts_cartesian[:, 2] - + pts_cartesian[:, 0] = r * np.cos(phi) + self.origin[0] + pts_cartesian[:, 1] = r * np.sin(phi) + self.origin[1] + pts_cartesian[:, 2] += self.origin[2] + return super().write_data_to_vtk( points=pts_cartesian, filename=filename, From c71c13496356a0e65cc1ac284608a31845f0f934 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Feb 2023 16:32:05 -0500 Subject: [PATCH 54/58] added grids in statepoint specs --- docs/source/io_formats/statepoint.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 9a08eed737..c0db5ea923 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -78,6 +78,12 @@ The current version of the statepoint file format is 17.0. of mesh. - **width** (*double[]*) -- Width of each mesh cell in each dimension. + - **Cylindrical & Spherical Mesh Only:** + - **r_grid** (*double[]*) -- The mesh divisions along the r-axis. + - **phi_grid** (*double[]*) -- The mesh divisions along the phi-axis. + - **origin** (*double[]*) -- The origin in cartesian coordinates. + - **Spherical Mesh Only:** + - **theta_grid** (*double[]*) -- The mesh divisions along the theta-axis. - **Unstructured Mesh Only:** - **filename** (*char[]*) -- Name of the mesh file. - **library** (*char[]*) -- Mesh library used to represent the From 59467633230660270fb74d0e6ceb0442b93e9680 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 4 Mar 2023 15:51:56 -0600 Subject: [PATCH 55/58] Moving origin data member down into Cylindrical/Spherical mesh using intermediate class --- include/openmc/mesh.h | 30 ++++++++++++++++++++++++++---- src/mesh.cpp | 14 ++++++++------ 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 9554a51436..6ae0761065 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -75,6 +75,12 @@ public: // Methods + //! Update a position to the local coordinates of the mesh + virtual void to_local_coords(Position& r) const {}; + + //! Return a position in the local coordinates of the mesh + virtual Position local_coords(const Position& r) const { return r; }; + //! Determine which bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -160,7 +166,7 @@ public: } }; - int get_bin(Position r) const override; + virtual int get_bin(Position r) const; int n_bins() const override; @@ -241,11 +247,27 @@ public: xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh std::array shape_; //!< Number of mesh elements in each dimension - Position origin_ {0.0, 0.0, 0.0}; protected: }; +class PeriodicStructuredMesh : public StructuredMesh { + +public: + PeriodicStructuredMesh() = default; + PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; + + void to_local_coords(Position& r) const override { r -= origin_; }; + + Position local_coords(const Position& r) const override + { + return r - origin_; + }; + + // Data members + Position origin_ {0.0, 0.0, 0.0}; +}; + //============================================================================== //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes //============================================================================== @@ -336,7 +358,7 @@ public: int set_grid(); }; -class CylindricalMesh : public StructuredMesh { +class CylindricalMesh : public PeriodicStructuredMesh { public: // Constructors CylindricalMesh() = default; @@ -390,7 +412,7 @@ private: } }; -class SphericalMesh : public StructuredMesh { +class SphericalMesh : public PeriodicStructuredMesh { public: // Constructors SphericalMesh() = default; diff --git a/src/mesh.cpp b/src/mesh.cpp index e0fa305edd..eca1d2a3fe 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -446,8 +446,8 @@ void StructuredMesh::raytrace_mesh( // translate start and end positions, // this needs to come after the get_indices call because it does its own translation - r0 -= origin_; - r1 -= origin_; + to_local_coords(r0); + to_local_coords(r1); // Calculate initial distances to next surfaces in all three dimensions std::array distances; @@ -952,7 +952,8 @@ void RectilinearMesh::to_hdf5(hid_t group) const // CylindricalMesh implementation //============================================================================== -CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} +CylindricalMesh::CylindricalMesh(pugi::xml_node node) + : PeriodicStructuredMesh {node} { n_dimension_ = 3; @@ -976,7 +977,7 @@ std::string CylindricalMesh::get_mesh_type() const StructuredMesh::MeshIndex CylindricalMesh::get_indices( Position r, bool& in_mesh) const { - r -= origin_; + to_local_coords(r); Position mapped_r; mapped_r[0] = std::hypot(r.x, r.y); @@ -1189,7 +1190,8 @@ void CylindricalMesh::to_hdf5(hid_t group) const // SphericalMesh implementation //============================================================================== -SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} +SphericalMesh::SphericalMesh(pugi::xml_node node) + : PeriodicStructuredMesh {node} { n_dimension_ = 3; @@ -1213,7 +1215,7 @@ std::string SphericalMesh::get_mesh_type() const StructuredMesh::MeshIndex SphericalMesh::get_indices( Position r, bool& in_mesh) const { - r -= origin_; + to_local_coords(r); Position mapped_r; mapped_r[0] = r.norm(); From b02e558fae56486aed8d97ee1f00e8831be3f317 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 4 Mar 2023 23:57:55 -0600 Subject: [PATCH 56/58] Adding docstring. Updating method name --- include/openmc/mesh.h | 6 +++--- src/mesh.cpp | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 6ae0761065..f19a22c09f 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -76,7 +76,7 @@ public: // Methods //! Update a position to the local coordinates of the mesh - virtual void to_local_coords(Position& r) const {}; + virtual void local_coords(Position& r) const {}; //! Return a position in the local coordinates of the mesh virtual Position local_coords(const Position& r) const { return r; }; @@ -257,7 +257,7 @@ public: PeriodicStructuredMesh() = default; PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; - void to_local_coords(Position& r) const override { r -= origin_; }; + void local_coords(Position& r) const override { r -= origin_; }; Position local_coords(const Position& r) const override { @@ -265,7 +265,7 @@ public: }; // Data members - Position origin_ {0.0, 0.0, 0.0}; + Position origin_ {0.0, 0.0, 0.0}; //!< Origin of the mesh }; //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index eca1d2a3fe..9020a95aaf 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -446,8 +446,8 @@ void StructuredMesh::raytrace_mesh( // translate start and end positions, // this needs to come after the get_indices call because it does its own translation - to_local_coords(r0); - to_local_coords(r1); + local_coords(r0); + local_coords(r1); // Calculate initial distances to next surfaces in all three dimensions std::array distances; @@ -977,7 +977,7 @@ std::string CylindricalMesh::get_mesh_type() const StructuredMesh::MeshIndex CylindricalMesh::get_indices( Position r, bool& in_mesh) const { - to_local_coords(r); + local_coords(r); Position mapped_r; mapped_r[0] = std::hypot(r.x, r.y); @@ -1215,7 +1215,7 @@ std::string SphericalMesh::get_mesh_type() const StructuredMesh::MeshIndex SphericalMesh::get_indices( Position r, bool& in_mesh) const { - to_local_coords(r); + local_coords(r); Position mapped_r; mapped_r[0] = r.norm(); From ab4237c037a3f142e50168690c2c0dca5edefcf5 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 7 Mar 2023 14:36:07 -0500 Subject: [PATCH 57/58] use get_reshaped_data in tests --- tests/unit_tests/test_cylindrical_mesh.py | 2 +- tests/unit_tests/test_spherical_mesh.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index e67f21b15c..9c8ce79f91 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -93,6 +93,6 @@ def test_offset_mesh(model, estimator, origin): # check that the lower half of the mesh contains a tally result # and that the upper half is zero - mean = tally.mean.reshape(mesh.dimension[::-1]).T + mean = tally.get_reshaped_data('mean', expand_dims=True) # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 # assert np.count_nonzero(mean[:, :, 5:]) == 0 diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index d6ac435c2e..489c4bd34e 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -93,6 +93,6 @@ def test_offset_mesh(model, estimator, origin): # check that the lower half of the mesh contains a tally result # and that the upper half is zero - mean = tally.mean.reshape(mesh.dimension[::-1]).T + mean = tally.get_reshaped_data('mean', expand_dims=True) # assert np.count_nonzero(mean[:, :, :5]) == mean.size / 2 # assert np.count_nonzero(mean[:, :, 5:]) == 0 From bb5f727f85325cc2d4e94d00b1940c9d744dd568 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 21 Mar 2023 13:48:22 -0500 Subject: [PATCH 58/58] Update IO format docs for meshes --- docs/source/io_formats/statepoint.rst | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index c0db5ea923..09b1b11442 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -72,12 +72,17 @@ The current version of the statepoint file format is 17.0. :Datasets: - **type** (*char[]*) -- Type of mesh. - **dimension** (*int*) -- Number of mesh cells in each dimension. - - **lower_left** (*double[]*) -- Coordinates of lower-left corner of - mesh. - - **upper_right** (*double[]*) -- Coordinates of upper-right corner - of mesh. - - **width** (*double[]*) -- Width of each mesh cell in each - dimension. + - **Regular Mesh Only:** + - **lower_left** (*double[]*) -- Coordinates of lower-left corner of + mesh. + - **upper_right** (*double[]*) -- Coordinates of upper-right corner + of mesh. + - **width** (*double[]*) -- Width of each mesh cell in each + dimension. + - **Rectilinear Mesh Only:** + - **x_grid** (*double[]*) -- Mesh divisions along the x-axis. + - **y_grid** (*double[]*) -- Mesh divisions along the y-axis. + - **z_grid** (*double[]*) -- Mesh divisions along the z-axis. - **Cylindrical & Spherical Mesh Only:** - **r_grid** (*double[]*) -- The mesh divisions along the r-axis. - **phi_grid** (*double[]*) -- The mesh divisions along the phi-axis.