From 687d43b0ab6e3e3ddfe867767c210f0a1cfb38a7 Mon Sep 17 00:00:00 2001 From: Patrick Myers Date: Wed, 25 May 2022 09:06:53 -0500 Subject: [PATCH 01/71] added external xtensor within PhotonInteraction for all shell photoelectric xs to utilize contiguous memory in calculate_xs --- include/openmc/photon.h | 2 +- src/photon.cpp | 26 +++++++++++--------------- src/physics.cpp | 6 +++--- 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index e8ec82d107..50ee228c44 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -36,7 +36,6 @@ public: int threshold; double n_electrons; double binding_energy; - xt::xtensor cross_section; vector transitions; }; @@ -82,6 +81,7 @@ public: // Photoionization and atomic relaxation data vector shells_; + xt::xtensor cross_sections_; // Compton profile data xt::xtensor profile_pdf_; diff --git a/src/photon.cpp b/src/photon.cpp index b0cd8aa87d..3f9ff6dbf3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -16,6 +16,9 @@ #include "xtensor/xbuilder.hpp" #include "xtensor/xoperation.hpp" #include "xtensor/xview.hpp" +#include "xtensor/xmath.hpp" +#include "xtensor/xslice.hpp" +#include "xtensor/xtensor_forward.hpp" #include #include @@ -125,6 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) } shells_.resize(n_shell); + cross_sections_.resize({energy_.size(), n_shell}); // Create mapping from designator to index std::unordered_map shell_map; @@ -155,13 +159,14 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_attribute(tgroup, "num_electrons", shell.n_electrons); // Read subshell cross section + xt::xtensor xs; dset = open_dataset(tgroup, "xs"); read_attribute(dset, "threshold_idx", shell.threshold); close_dataset(dset); - read_dataset(tgroup, "xs", shell.cross_section); + read_dataset(tgroup, "xs", xs); - auto& xs = shell.cross_section; - xs = xt::where(xs > 0.0, xt::log(xs), -500.0); + auto cross_section = xt::view(cross_sections_, xt::range(shell.threshold, shell.threshold + xs.size()), i); + cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); if (object_exists(tgroup, "transitions")) { // Determine dimensions of transitions @@ -565,19 +570,10 @@ void PhotonInteraction::calculate_xs(Particle& p) const incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid))); // Calculate microscopic photoelectric cross section - xs.photoelectric = 0.0; - for (const auto& shell : shells_) { - // Check threshold of reaction - int i_start = shell.threshold; - if (i_grid < i_start) - continue; + const auto& xs_upper = xt::row(cross_sections_, i_grid); + const auto& xs_lower = xt::row(cross_sections_, i_grid + 1); - // Evaluation subshell photoionization cross section - xs.photoelectric += - std::exp(shell.cross_section(i_grid - i_start) + - f * (shell.cross_section(i_grid + 1 - i_start) - - shell.cross_section(i_grid - i_start))); - } + xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0]; // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 13dd063de5..38d74dbede 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -358,9 +358,9 @@ void sample_photon_reaction(Particle& p) continue; // Evaluation subshell photoionization cross section - double xs = std::exp(shell.cross_section(i_grid - i_start) + - f * (shell.cross_section(i_grid + 1 - i_start) - - shell.cross_section(i_grid - i_start))); + double xs = std::exp(element.cross_sections_(i_grid - i_start) + + f * (element.cross_sections_(i_grid + 1 - i_start) - + element.cross_sections_(i_grid - i_start))); prob += xs; if (prob > cutoff) { From 6df0286b2e20e39b4eff8d93b262821086a9304f Mon Sep 17 00:00:00 2001 From: Patrick Myers Date: Wed, 25 May 2022 09:14:10 -0500 Subject: [PATCH 02/71] added clang-format --- src/photon.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 3f9ff6dbf3..1be515a621 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -14,11 +14,11 @@ #include "openmc/settings.h" #include "xtensor/xbuilder.hpp" -#include "xtensor/xoperation.hpp" -#include "xtensor/xview.hpp" #include "xtensor/xmath.hpp" +#include "xtensor/xoperation.hpp" #include "xtensor/xslice.hpp" #include "xtensor/xtensor_forward.hpp" +#include "xtensor/xview.hpp" #include #include @@ -165,7 +165,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view(cross_sections_, xt::range(shell.threshold, shell.threshold + xs.size()), i); + auto cross_section = xt::view(cross_sections_, + xt::range(shell.threshold, shell.threshold + xs.size()), i); cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); if (object_exists(tgroup, "transitions")) { From 61db44f40e33d132d895974852c448beeacc6770 Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 9 Jun 2022 09:52:50 -0500 Subject: [PATCH 03/71] Xtensor for photoelectric cross sections matches the flux spectrum for uranium of the previous iteration of OpenMC --- src/photon.cpp | 15 +++++++++------ src/physics.cpp | 22 +++++++++++----------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 1be515a621..bb715cdeaf 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -128,7 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) } shells_.resize(n_shell); - cross_sections_.resize({energy_.size(), n_shell}); + cross_sections_ = xt::zeros({energy_.size(), n_shell}); // Create mapping from designator to index std::unordered_map shell_map; @@ -166,8 +166,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_dataset(tgroup, "xs", xs); auto cross_section = xt::view(cross_sections_, - xt::range(shell.threshold, shell.threshold + xs.size()), i); - cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); + xt::range(shell.threshold, xt::placeholders::_), i); + cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { // Determine dimensions of transitions @@ -571,10 +571,13 @@ void PhotonInteraction::calculate_xs(Particle& p) const incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid))); // Calculate microscopic photoelectric cross section - const auto& xs_upper = xt::row(cross_sections_, i_grid); - const auto& xs_lower = xt::row(cross_sections_, i_grid + 1); + xs.photoelectric = 0.0; + const auto& xs_lower = xt::row(cross_sections_, i_grid); + const auto& xs_upper = xt::row(cross_sections_, i_grid + 1); - xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0]; + for (int i = 0; i < xs_upper.size(); ++i) + if (xs_lower(i) != 0) + xs.photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 38d74dbede..26b3a3520a 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -29,6 +29,7 @@ #include // for max, min, max_element #include // for sqrt, exp, log, abs, copysign +#include namespace openmc { @@ -344,25 +345,24 @@ void sample_photon_reaction(Particle& p) // Photoelectric effect double prob_after = prob + micro.photoelectric; + int i_grid = micro.index_grid; + double f = micro.interp_factor; + const auto& xs_lower = xt::row(element.cross_sections_, i_grid); + const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); + if (prob_after > cutoff) { for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { const auto& shell {element.shells_[i_shell]}; - // Get grid index and interpolation factor - int i_grid = micro.index_grid; - double f = micro.interp_factor; - // Check threshold of reaction - int i_start = shell.threshold; - if (i_grid < i_start) + if (xs_lower(i_shell) == 0) continue; - // Evaluation subshell photoionization cross section - double xs = std::exp(element.cross_sections_(i_grid - i_start) + - f * (element.cross_sections_(i_grid + 1 - i_start) - - element.cross_sections_(i_grid - i_start))); + // Evaluation subshell photoionization cross section + prob += std::exp(xs_lower(i_shell) + + f * (xs_upper(i_shell) - + xs_lower(i_shell))); - prob += xs; if (prob > cutoff) { double E_electron = p.E() - shell.binding_energy; From 70b1706460166ff6d0a05744ad804ef9fd0671ef Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 9 Jun 2022 10:39:46 -0500 Subject: [PATCH 04/71] Moved initialized variables into if statement --- src/physics.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 26b3a3520a..0df8e63c0a 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -345,12 +345,14 @@ void sample_photon_reaction(Particle& p) // Photoelectric effect double prob_after = prob + micro.photoelectric; - int i_grid = micro.index_grid; - double f = micro.interp_factor; - const auto& xs_lower = xt::row(element.cross_sections_, i_grid); - const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); if (prob_after > cutoff) { + + int i_grid = micro.index_grid; + double f = micro.interp_factor; + const auto& xs_lower = xt::row(element.cross_sections_, i_grid); + const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); + for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { const auto& shell {element.shells_[i_shell]}; From 98a7d5eabc972a986a64307bc709dda7ecadaf3d Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 10 Jun 2022 11:57:35 -0500 Subject: [PATCH 05/71] added clang format --- src/photon.cpp | 7 ++++--- src/physics.cpp | 5 ++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index bb715cdeaf..85bd435c99 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -165,8 +165,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view(cross_sections_, - xt::range(shell.threshold, xt::placeholders::_), i); + auto cross_section = xt::view( + cross_sections_, xt::range(shell.threshold, xt::placeholders::_), i); cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { @@ -577,7 +577,8 @@ void PhotonInteraction::calculate_xs(Particle& p) const for (int i = 0; i < xs_upper.size(); ++i) if (xs_lower(i) != 0) - xs.photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); + xs.photoelectric += + std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 0df8e63c0a..cc712affd7 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -361,9 +361,8 @@ void sample_photon_reaction(Particle& p) continue; // Evaluation subshell photoionization cross section - prob += std::exp(xs_lower(i_shell) + - f * (xs_upper(i_shell) - - xs_lower(i_shell))); + prob += std::exp( + xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); if (prob > cutoff) { double E_electron = p.E() - shell.binding_energy; From a4f98b489c161ed4b7ac26cca4e945d535ccaba3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 23 May 2022 14:20:01 -0500 Subject: [PATCH 06/71] Change values in dict returned by Material.get_nuclide_atom_densities --- openmc/cell.py | 8 ++++---- openmc/deplete/operator.py | 15 +++++++-------- openmc/material.py | 20 ++++++++++++-------- openmc/plotter.py | 12 ++++-------- tests/unit_tests/test_deplete_activation.py | 2 +- tests/unit_tests/test_material.py | 15 +++++++-------- tests/unit_tests/test_model.py | 18 ++++++++---------- 7 files changed, 43 insertions(+), 47 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 103a34fe63..4b419e1c6e 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -214,8 +214,8 @@ class Cell(IDManagerMixin): self._atoms = self._fill.get_nuclide_atom_densities() # Convert to total number of atoms - for key, nuclide in self._atoms.items(): - atom = nuclide[1] * self._volume * 1.0e+24 + for key, atom_per_bcm in self._atoms.items(): + atom = atom_per_bcm * self._volume * 1.0e+24 self._atoms[key] = atom elif self.fill_type == 'distribmat': @@ -224,12 +224,12 @@ class Cell(IDManagerMixin): partial_volume = self.volume / len(self.fill) self._atoms = OrderedDict() for mat in self.fill: - for key, nuclide in mat.get_nuclide_atom_densities().items(): + for key, atom_per_bcm in mat.get_nuclide_atom_densities().items(): # To account for overlap of nuclides between distribmat # we need to append new atoms to any existing value # hence it is necessary to ask for default. atom = self._atoms.setdefault(key, 0) - atom += nuclide[1] * partial_volume * 1.0e+24 + atom += atom_per_bcm * partial_volume * 1.0e+24 self._atoms[key] = atom else: diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 327c762e08..7ef874e87f 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -528,9 +528,9 @@ class Operator(TransportOperator): """ mat_id = str(mat.id) - for nuclide, density in mat.get_nuclide_atom_densities().values(): - number = density * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, number) + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) def _set_number_from_results(self, mat, prev_res): """Extracts material nuclides and number densities. @@ -556,16 +556,15 @@ class Operator(TransportOperator): # Merge lists of nuclides, with the same order for every calculation geom_nuc_densities.update(depl_nuc) - for nuclide in geom_nuc_densities.keys(): + for nuclide, atom_per_bcm in geom_nuc_densities.items(): if nuclide in depl_nuc: concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] volume = prev_res[-1].volume[mat_id] - number = concentration / volume + atom_per_cc = concentration / volume else: - density = geom_nuc_densities[nuclide][1] - number = density * 1.0e24 + atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, number) + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) def initial_condition(self): """Performs final setup and returns initial condition. diff --git a/openmc/material.py b/openmc/material.py index 3c2a98ffd3..073dab7121 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -259,10 +259,10 @@ class Material(IDManagerMixin): if self.volume is None: raise ValueError("Volume must be set in order to determine mass.") density = 0.0 - for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items(): Z = openmc.data.zam(nuc)[0] if Z >= 90: - density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ + density += 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \ / openmc.data.AVOGADRO return density*self.volume @@ -786,11 +786,15 @@ class Material(IDManagerMixin): """Returns all nuclides in the material and their atomic densities in units of atom/b-cm + .. versionchanged:: 0.13.1 + The values in the dictionary were changed from a tuple containing + the nuclide name and the density to just the density. + Returns ------- nuclides : dict - Dictionary whose keys are nuclide names and values are tuples of - (nuclide, density in atom/b-cm) + Dictionary whose keys are nuclide names and values are densities in + [atom/b-cm] """ @@ -850,7 +854,7 @@ class Material(IDManagerMixin): nuclides = OrderedDict() for n, nuc in enumerate(nucs): - nuclides[nuc] = (nuc, nuc_densities[n]) + nuclides[nuc] = nuc_densities[n] return nuclides @@ -870,9 +874,9 @@ class Material(IDManagerMixin): """ mass_density = 0.0 - for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items(): if nuclide is None or nuclide == nuc: - density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ + density_i = 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \ / openmc.data.AVOGADRO mass_density += density_i return mass_density @@ -1092,7 +1096,7 @@ class Material(IDManagerMixin): nuclides_per_cc = defaultdict(float) mass_per_cc = defaultdict(float) for mat, wgt in zip(materials, wgts): - for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().items(): nuc_per_cc = wgt*1.e24*atoms_per_bcm nuclides_per_cc[nuc] += nuc_per_cc mass_per_cc[nuc] += nuc_per_cc*openmc.data.atomic_mass(nuc) / \ diff --git a/openmc/plotter.py b/openmc/plotter.py index 92acaf5bb5..16760868e3 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -542,15 +542,11 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., if isinstance(this, openmc.Material): # Expand elements in to nuclides with atomic densities - nuclides = this.get_nuclide_atom_densities() - # For ease of processing split out the nuclide and its fraction - nuc_fractions = {nuclide[1][0]: nuclide[1][1] - for nuclide in nuclides.items()} + nuc_fractions = this.get_nuclide_atom_densities() # Create a dict of [nuclide name] = nuclide object to carry forward # with a common nuclides format between openmc.Material and # openmc.Element objects - nuclides = {nuclide[1][0]: nuclide[1][0] - for nuclide in nuclides.items()} + nuclides = {nuclide: nuclide for nuclide in nuc_fractions} else: # Expand elements in to nuclides with atomic densities nuclides = this.expand(1., 'ao', enrichment=enrichment, @@ -885,13 +881,13 @@ def _calculate_mgxs_elem_mat(this, types, library, orders=None, # Check to see if we have nuclides/elements or a macroscopic object if this._macroscopic is not None: # We have macroscopics - nuclides = {this._macroscopic: (this._macroscopic, this.density)} + nuclides = {this._macroscopic: this.density} else: # Expand elements in to nuclides with atomic densities nuclides = this.get_nuclide_atom_densities() # For ease of processing split out nuc and nuc_density - nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()] + nuc_fraction = list(nuclides.values()) else: T = temperature # Expand elements in to nuclides with atomic densities diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 0b82a5fbcf..cb3b86b9c8 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -92,7 +92,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts w = model.geometry.get_materials_by_name('tungsten')[0] atom_densities = w.get_nuclide_atom_densities() - atom_per_cc = 1e24 * atom_densities['W186'][1] # Density in atom/cm^3 + atom_per_cc = 1e24 * atom_densities['W186'] # Density in atom/cm^3 n0 = atom_per_cc * w.volume # Absolute number of atoms # Pick a random irradiation time and then determine necessary source rate to diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index c15a4b9e55..374b5a3d8a 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -100,7 +100,7 @@ def test_add_elements_by_formula(): m.add_elements_from_formula('Li4SiO4') # checking the ratio of elements is 4:1:4 for Li:Si:O elem = defaultdict(float) - for nuclide, adens in m.get_nuclide_atom_densities().values(): + for nuclide, adens in m.get_nuclide_atom_densities().items(): if nuclide.startswith("Li"): elem["Li"] += adens if nuclide.startswith("Si"): @@ -117,7 +117,7 @@ def test_add_elements_by_formula(): 'O16': 0.443386, 'O17': 0.000168} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing the correct nuclides are added to the Material when enriched m = openmc.Material() @@ -129,7 +129,7 @@ def test_add_elements_by_formula(): 'O16': 0.443386, 'O17': 0.000168} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing the use of brackets m = openmc.Material() @@ -137,7 +137,7 @@ def test_add_elements_by_formula(): # checking the ratio of elements is 2:2:6 for Mg:N:O elem = defaultdict(float) - for nuclide, adens in m.get_nuclide_atom_densities().values(): + for nuclide, adens in m.get_nuclide_atom_densities().items(): if nuclide.startswith("Mg"): elem["Mg"] += adens if nuclide.startswith("N"): @@ -155,7 +155,7 @@ def test_add_elements_by_formula(): 'O16': 0.599772, 'O17': 0.000227} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing non integer multiplier results in a value error m = openmc.Material() @@ -289,8 +289,7 @@ def test_get_nuclide_densities(uo2): def test_get_nuclide_atom_densities(uo2): - nucs = uo2.get_nuclide_atom_densities() - for nuc, density in nucs.values(): + for nuc, density in uo2.get_nuclide_atom_densities().items(): assert nuc in ('U235', 'O16') assert density > 0 @@ -341,7 +340,7 @@ def test_borated_water(): 'O16':2.4672e-02} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) assert m.id == 50 # Test the Celsius conversion. diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 10417b20ed..9b05ed2a5c 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -391,16 +391,14 @@ def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm): assert openmc.lib.materials[1].get_density('atom/b-cm') == \ pytest.approx(0.06891296988603757, abs=1e-13) mat_a_dens = np.sum( - [v[1] for v in test_model.materials[0]. - get_nuclide_atom_densities().values()]) + list(test_model.materials[0].get_nuclide_atom_densities().values())) assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8) # Change the density test_model.update_densities(['UO2'], 2.) assert openmc.lib.materials[1].get_density('atom/b-cm') == \ pytest.approx(2., abs=1e-13) mat_a_dens = np.sum( - [v[1] for v in test_model.materials[0]. - get_nuclide_atom_densities().values()]) + list(test_model.materials[0].get_nuclide_atom_densities().values())) assert mat_a_dens == pytest.approx(2., abs=1e-8) # Now lets do the cell temperature updates. @@ -441,7 +439,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): test_model = openmc.Model(geom, mats, settings, tals, plots) initial_mat = mats[0].clone() - initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1] + initial_u = initial_mat.get_nuclide_atom_densities()['U235'] # Note that the chain file includes only U-235 fission to a stable Xe136 w/ # a yield of 100%. Thus all the U235 we lose becomes Xe136 @@ -453,8 +451,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): operator_kwargs=op_kwargs, power=1., output=False) # Get the new Xe136 and U235 atom densities - after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] - after_u = mats[0].get_nuclide_atom_densities()['U235'][1] + after_xe = mats[0].get_nuclide_atom_densities()['Xe136'] + after_u = mats[0].get_nuclide_atom_densities()['U235'] assert after_xe + after_u == pytest.approx(initial_u, abs=1e-15) assert test_model.is_initialized is False @@ -462,7 +460,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): mats[0].nuclides.clear() densities = initial_mat.get_nuclide_atom_densities() tot_density = 0. - for nuc, density in densities.values(): + for nuc, density in densities.items(): mats[0].add_nuclide(nuc, density) tot_density += density mats[0].set_density('atom/b-cm', tot_density) @@ -473,8 +471,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): operator_kwargs=op_kwargs, power=1., output=False) # Get the new Xe136 and U235 atom densities - after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] - after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1] + after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'] + after_lib_u = mats[0].get_nuclide_atom_densities()['U235'] assert after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15) assert test_model.is_initialized is True From 7fcc20472690c49cf2b0508b099dc440b6249dc4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 23 May 2022 14:43:34 -0500 Subject: [PATCH 07/71] Implement Material.get_nuclide_atoms method --- openmc/material.py | 19 +++++++++++++++++++ tests/unit_tests/test_material.py | 10 ++++++++++ 2 files changed, 29 insertions(+) diff --git a/openmc/material.py b/openmc/material.py index 073dab7121..ea42d64b16 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -858,6 +858,25 @@ class Material(IDManagerMixin): return nuclides + def get_nuclide_atoms(self): + """Return number of atoms of each nuclide in the material + + .. versionadded:: 0.13.1 + + Returns + ------- + dict + Dictionary whose keys are nuclide names and values are number of + atoms present in the material. + + """ + if self.volume is None: + raise ValueError("Volume must be set in order to determine atoms.") + atoms = {} + for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items(): + atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume + return atoms + def get_mass_density(self, nuclide=None): """Return mass density of one or all nuclides diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 374b5a3d8a..fd0e4f6247 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -294,6 +294,16 @@ def test_get_nuclide_atom_densities(uo2): assert density > 0 +def test_get_nuclide_atoms(): + mat = openmc.Material() + mat.add_nuclide('Li6', 1.0) + mat.set_density('atom/cm3', 3.26e20) + mat.volume = 100.0 + + atoms = mat.get_nuclide_atoms() + assert atoms['Li6'] == pytest.approx(mat.density * mat.volume) + + def test_mass(): m = openmc.Material() m.add_nuclide('Zr90', 1.0, 'wo') From 19494cccd022cdf9ee9c0c8e6461bea494d95e69 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 24 May 2022 12:47:36 -0500 Subject: [PATCH 08/71] Add check in Decay.decay_constant for erroneous half-life --- openmc/data/decay.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 5ca2b235a5..51a594bc80 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -465,11 +465,10 @@ class Decay(EqualityMixin): @property def decay_constant(self): - if hasattr(self.half_life, 'n'): - return log(2.)/self.half_life - else: - mu, sigma = self.half_life - return ufloat(log(2.)/mu, log(2.)/mu**2*sigma) + if self.half_life.n == 0.0: + name = self.nuclide['name'] + raise ValueError(f"{name} is listed as unstable but has a zero half-life.") + return log(2.)/self.half_life @property def decay_energy(self): From 3c20b1e415e89b1b2b2b8e4953710f0b6851fb62 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Jun 2022 08:00:16 -0500 Subject: [PATCH 09/71] Fix use of get_nuclide_atom_densities in Material.activity property --- openmc/material.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index ea42d64b16..1341746dc0 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -86,6 +86,9 @@ class Material(IDManagerMixin): fissionable_mass : float Mass of fissionable nuclides in the material in [g]. Requires that the :attr:`volume` attribute is set. + activity : float + Activity of the material in [Bq]. Requires that the :attr:`volume` + attribute is set. """ @@ -152,7 +155,7 @@ class Material(IDManagerMixin): for key, value in atoms_per_barn_cm.items(): half_life = openmc.data.half_life(key) if half_life: - total_activity += value[1] / half_life + total_activity += value / half_life total_activity *= math.log(2) * 1e24 * self.volume return total_activity From fc615451d04c358d0f32bb72837fc52992fb39b0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 15:54:49 +0200 Subject: [PATCH 10/71] added method to Tally --- openmc/tallies.py | 86 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index d0355f14ee..f559c392cc 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -8,6 +8,8 @@ from pathlib import Path from xml.etree import ElementTree as ET import h5py +import vtk +import vtk.util.numpy_support as nps import numpy as np import pandas as pd import scipy.sparse as sps @@ -457,6 +459,54 @@ class Tally(IDManagerMixin): self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) self._sparse = False + def write_to_vtk(self, filename): + mesh = None + for f in self.filters: + if isinstance(f, openmc.MeshFilter): + if isinstance(f.mesh, (openmc.RegularMesh, openmc.CylindricalMesh)): + mesh = f.mesh + break + + if not mesh: + raise ValueError( + "write_to_vtk only works with openmc.RegularMesh, openmc.CylindricalMesh" + ) + + if isinstance(mesh, openmc.RegularMesh): + vtk_grid = voxels_to_vtk( + x_vals=np.linspace( + mesh.lower_left[0], + mesh.upper_right[0], + num=mesh.dimension[0] + 1, + ), + y_vals=np.linspace( + mesh.lower_left[1], + mesh.upper_right[1], + num=mesh.dimension[1] + 1, + ), + z_vals=np.linspace( + mesh.lower_left[2], + mesh.upper_right[2], + num=mesh.dimension[2] + 1, + ), + mean=self.mean, + std_dev=self.std_dev, + cylindrical=False, + ) + elif isinstance(mesh, openmc.CylindricalMesh): + vtk_grid = voxels_to_vtk( + x_vals=mesh.r_grid, + y_vals=mesh.phi_grid, + z_vals=mesh.z_grid, + mean=self.mean, + std_dev=self.std_dev, + cylindrical=True, + ) + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + def remove_score(self, score): """Remove a score from the tally @@ -3229,3 +3279,39 @@ class Tallies(cv.CheckedList): tallies.append(tally) return cls(tallies) + + +def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): + vtk_box = vtk.vtkStructuredGrid() + + vtk_box.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # points + points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + if cylindrical: + points_cartesian = np.copy(points) + r = points[:, 0] + phi = points[:, 1] + z = points[:, 2] + points_cartesian[:, 0] = r * np.cos(phi) + points_cartesian[:, 1] = r * np.sin(phi) + points_cartesian[:, 2] = z + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + + vtk_box.SetPoints(vtkPts) + + # # data + mean_array = vtk.vtkDoubleArray() + mean_array.SetName("mean") + mean_array.SetArray(mean, mean.size, True) + + std_dev_array = vtk.vtkDoubleArray() + std_dev_array.SetName("std_dev") + std_dev_array.SetArray(std_dev, std_dev.size, True) + + vtk_box.GetCellData().AddArray(mean_array) + vtk_box.GetCellData().AddArray(std_dev_array) + + return vtk_box From e870ab2040e6d2a9894eb99e413a5174f83ae546 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:00:54 +0200 Subject: [PATCH 11/71] added RectilinearMesh --- openmc/tallies.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index f559c392cc..20e6e39964 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -463,13 +463,13 @@ class Tally(IDManagerMixin): mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): - if isinstance(f.mesh, (openmc.RegularMesh, openmc.CylindricalMesh)): + if isinstance(f.mesh, (openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh)): mesh = f.mesh break if not mesh: raise ValueError( - "write_to_vtk only works with openmc.RegularMesh, openmc.CylindricalMesh" + "write_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" ) if isinstance(mesh, openmc.RegularMesh): @@ -493,6 +493,15 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=False, ) + if isinstance(mesh, openmc.RectilinearMesh): + vtk_grid = voxels_to_vtk( + x_vals=mesh.x_grid, + y_vals=mesh.y_grid, + z_vals=mesh.z_grid, + mean=self.mean, + std_dev=self.std_dev, + cylindrical=False, + ) elif isinstance(mesh, openmc.CylindricalMesh): vtk_grid = voxels_to_vtk( x_vals=mesh.r_grid, From 06c5d3ec584c9ef9a62b08c561945b34e2e6b4e9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:04:33 +0200 Subject: [PATCH 12/71] added comments --- openmc/tallies.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 20e6e39964..6797253e1c 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -460,6 +460,7 @@ class Tally(IDManagerMixin): self._sparse = False def write_to_vtk(self, filename): + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): @@ -511,6 +512,8 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=True, ) + + # write the .vtk file writer = vtk.vtkStructuredGridWriter() writer.SetFileName(filename) writer.SetInputData(vtk_grid) @@ -3291,36 +3294,32 @@ class Tallies(cv.CheckedList): def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): - vtk_box = vtk.vtkStructuredGrid() + vtk_grid = vtk.vtkStructuredGrid() - vtk_box.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - # points + # create points points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if cylindrical: + if cylindrical: # transform points to cartesian coordinates points_cartesian = np.copy(points) - r = points[:, 0] - phi = points[:, 1] - z = points[:, 2] + r, phi, z = points[:, 0], points[:, 1], points[:, 2] points_cartesian[:, 0] = r * np.cos(phi) points_cartesian[:, 1] = r * np.sin(phi) points_cartesian[:, 2] = z vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + vtk_grid.SetPoints(vtkPts) - vtk_box.SetPoints(vtkPts) - - # # data + # add mean and std dev data mean_array = vtk.vtkDoubleArray() mean_array.SetName("mean") mean_array.SetArray(mean, mean.size, True) + vtk_grid.GetCellData().AddArray(mean_array) std_dev_array = vtk.vtkDoubleArray() std_dev_array.SetName("std_dev") std_dev_array.SetArray(std_dev, std_dev.size, True) + vtk_grid.GetCellData().AddArray(std_dev_array) - vtk_box.GetCellData().AddArray(mean_array) - vtk_box.GetCellData().AddArray(std_dev_array) - - return vtk_box + return vtk_grid From 30b06fa2910a51ac04a249170bff27d54cb3f26a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:08:36 +0200 Subject: [PATCH 13/71] docstrings --- openmc/tallies.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 6797253e1c..77f0b8cd39 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -460,6 +460,15 @@ class Tally(IDManagerMixin): self._sparse = False def write_to_vtk(self, filename): + """Writes the tally to a vtk file + + Args: + filename (str): the filename (must end with .vtk) + + Raises: + ValueError: if no MeshFilter with appropriate mesh was found + (SphericalMesh not supported) + """ # check that tally has a MeshFilter mesh = None for f in self.filters: @@ -512,7 +521,7 @@ class Tally(IDManagerMixin): std_dev=self.std_dev, cylindrical=True, ) - + # write the .vtk file writer = vtk.vtkStructuredGridWriter() writer.SetFileName(filename) @@ -3294,6 +3303,20 @@ class Tallies(cv.CheckedList): def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): + """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. + + Args: + x_vals (list): X values. + y_vals (list): Y values. + z_vals (list): Z values. + mean (np.array): the tally mean. + std_dev (np.array): the tally standard deviation. + cylindrical (bool, optional): If set to True, cylindrical coordinates + (r, phi, z) are assumed. Defaults to True. + + Returns: + vtkStructuredGrid: a vtk object containing tally data on the appropriate grid + """ vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) From 06c6bcd7d269aa3b3317c8d64f9cf4318170468c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:09:32 +0200 Subject: [PATCH 14/71] added TODO --- openmc/tallies.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index 77f0b8cd39..c0f410524f 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3317,6 +3317,9 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid """ + + # TODO: should this be a method of Tally? + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) From 877a62c0af072e14c3b62991f21a6648352dc8bb Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:43:42 +0200 Subject: [PATCH 15/71] fixed cylindrical --- openmc/tallies.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index c0f410524f..2e31e8ee26 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3332,6 +3332,7 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): points_cartesian[:, 0] = r * np.cos(phi) points_cartesian[:, 1] = r * np.sin(phi) points_cartesian[:, 2] = z + points = points_cartesian vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) From 1eb8670f2b93cbbb42912dd7903c1d1ae18b1d33 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Jun 2022 16:46:10 +0200 Subject: [PATCH 16/71] mean and std_dev are np.zeros if None (for testing without running a sim) --- openmc/tallies.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/openmc/tallies.py b/openmc/tallies.py index 2e31e8ee26..75c7ea8f3b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3341,11 +3341,23 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): # add mean and std dev data mean_array = vtk.vtkDoubleArray() mean_array.SetName("mean") + if mean is None: + mean = np.zeros( + (len(x_vals) - 1) + * (len(y_vals) - 1) + * (len(z_vals) - 1) + ) mean_array.SetArray(mean, mean.size, True) vtk_grid.GetCellData().AddArray(mean_array) std_dev_array = vtk.vtkDoubleArray() std_dev_array.SetName("std_dev") + if std_dev is None: + std_dev = np.zeros( + (len(x_vals) - 1) + * (len(y_vals) - 1) + * (len(z_vals) - 1) + ) std_dev_array.SetArray(std_dev, std_dev.size, True) vtk_grid.GetCellData().AddArray(std_dev_array) From 3d0b404b7a99388296dda50d087d267bf3c6f917 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 28 Jun 2022 16:06:17 +0100 Subject: [PATCH 17/71] added multi stages option --- Dockerfile | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Dockerfile b/Dockerfile index 755d9f2a06..152590161a 100644 --- a/Dockerfile +++ b/Dockerfile @@ -15,7 +15,7 @@ # sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number -FROM debian:bullseye-slim +FROM debian:bullseye-slim AS dependencies # By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support ARG build_dagmc=off @@ -172,6 +172,8 @@ RUN if [ "$build_libmesh" = "on" ]; then \ && rm -rf ${LIBMESH_INSTALL_DIR}/build ${LIBMESH_INSTALL_DIR}/libmesh ; \ fi +FROM dependencies as build + # clone and install openmc RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \ @@ -207,5 +209,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && cd ../openmc && pip install .[test,depletion-mpi] \ && python -c "import openmc" +FROM build as release + # Download cross sections (NNDC and WMP) and ENDF data needed by test suite RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh From 419351bb260f837386bcd5c6e798405638820e17 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 28 Jun 2022 16:14:52 +0100 Subject: [PATCH 18/71] uppercase AS --- Dockerfile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index 152590161a..341390c9ae 100644 --- a/Dockerfile +++ b/Dockerfile @@ -172,7 +172,7 @@ RUN if [ "$build_libmesh" = "on" ]; then \ && rm -rf ${LIBMESH_INSTALL_DIR}/build ${LIBMESH_INSTALL_DIR}/libmesh ; \ fi -FROM dependencies as build +FROM dependencies AS build # clone and install openmc RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ @@ -209,7 +209,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && cd ../openmc && pip install .[test,depletion-mpi] \ && python -c "import openmc" -FROM build as release +FROM build AS release # Download cross sections (NNDC and WMP) and ENDF data needed by test suite RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh From 782d987fc836494166ffb66ab9e9f81888666582 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Tue, 28 Jun 2022 20:18:52 +0200 Subject: [PATCH 19/71] Optional vtk import Co-authored-by: Jonathan Shimwell --- openmc/tallies.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 75c7ea8f3b..b795fc83a2 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -8,8 +8,6 @@ from pathlib import Path from xml.etree import ElementTree as ET import h5py -import vtk -import vtk.util.numpy_support as nps import numpy as np import pandas as pd import scipy.sparse as sps @@ -470,6 +468,12 @@ class Tally(IDManagerMixin): (SphericalMesh not supported) """ # check that tally has a MeshFilter + try: + import vtk + import vtk.util.numpy_support as nps + except ImportError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ImportError(msg) mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From 7723ea492c12e17c2a8e2bef8f301386410b15a3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Jun 2022 09:45:38 -0500 Subject: [PATCH 20/71] Add Material.get_nuclide_activity method and decay_constant function --- docs/source/pythonapi/data.rst | 1 + openmc/data/data.py | 35 ++++++++++++++++++++++++++++-- openmc/material.py | 30 +++++++++++++++---------- tests/unit_tests/test_data_misc.py | 13 ++++++++++- 4 files changed, 64 insertions(+), 15 deletions(-) diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 6f5c006c24..287738774c 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -61,6 +61,7 @@ Core Functions atomic_mass atomic_weight + decay_constant dose_coefficients gnd_name half_life diff --git a/openmc/data/data.py b/openmc/data/data.py index 3a3cf58268..71b93293e7 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -3,7 +3,7 @@ import json import os import re from pathlib import Path -from math import sqrt +from math import sqrt, log from warnings import warn # Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions @@ -202,7 +202,7 @@ _GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)') # Used in half_life function as a cache _HALF_LIFE = {} - +_LOG_TWO = log(2.0) def atomic_mass(isotope): """Return atomic mass of isotope in atomic mass units. @@ -283,6 +283,8 @@ def half_life(isotope): Half-life values are from the `ENDF/B-VIII.0 decay sublibrary `_. + .. versionadded:: 0.13.1 + Parameters ---------- isotope : str @@ -302,6 +304,35 @@ def half_life(isotope): return _HALF_LIFE.get(isotope.lower()) + +def decay_constant(isotope): + """Return decay constant of isotope in [s^-1] + + Decay constants are based on half-life values from the + :func:`~openmc.data.half_life` function. When the isotope is stable, a decay + constant of zero is returned. + + .. versionadded:: 0.13.1 + + Parameters + ---------- + isotope : str + Name of isotope, e.g., 'Pu239' + + Returns + ------- + float + Decay constant of isotope in [s^-1] + + See also + -------- + openmc.data.half_life + + """ + t = half_life(isotope) + return _LOG_TWO / t if t else 0.0 + + def water_density(temperature, pressure=0.1013): """Return the density of liquid water at a given temperature and pressure. diff --git a/openmc/material.py b/openmc/material.py index 1341746dc0..e8285de7c2 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -3,7 +3,6 @@ from collections.abc import Iterable from copy import deepcopy from numbers import Real from pathlib import Path -import math import re import warnings from xml.etree import ElementTree as ET @@ -145,20 +144,10 @@ class Material(IDManagerMixin): return string - @property def activity(self): """Returns the total activity of the material in Becquerels.""" - - atoms_per_barn_cm = self.get_nuclide_atom_densities() - total_activity = 0 - for key, value in atoms_per_barn_cm.items(): - half_life = openmc.data.half_life(key) - if half_life: - total_activity += value / half_life - total_activity *= math.log(2) * 1e24 * self.volume - - return total_activity + return sum(self.get_nuclide_activity().values()) @property def name(self): @@ -861,6 +850,23 @@ class Material(IDManagerMixin): return nuclides + def get_nuclide_activity(self): + """Return activity in [Bq] for each nuclide in the material + + .. versionadded:: 0.13.1 + + Returns + ------- + dict + Dictionary whose keys are nuclide names and values are activity in + [Bq]. + """ + activity = {} + for nuclide, atoms in self.get_nuclide_atoms().items(): + inv_seconds = openmc.data.decay_constant(nuclide) + activity[nuclide] = inv_seconds * atoms + return activity + def get_nuclide_atoms(self): """Return number of atoms of each nuclide in the material diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index 6a81fb1e08..c3f7e3cfff 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -from collections.abc import Mapping +from math import log import os from pathlib import Path @@ -126,3 +126,14 @@ def test_zam(): assert openmc.data.zam('Am242_m10') == (95, 242, 10) with pytest.raises(ValueError): openmc.data.zam('garbage') + + +def test_half_life(): + assert openmc.data.half_life('H2') is None + assert openmc.data.half_life('U235') == pytest.approx(2.22102e16) + assert openmc.data.half_life('Am242') == pytest.approx(57672.0) + assert openmc.data.half_life('Am242_m1') == pytest.approx(4449622000.0) + assert openmc.data.decay_constant('H2') == 0.0 + assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16) + assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0) + assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0) From 6eaab920eeaddac452ad07e25015357c3658f139 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:18:10 +0200 Subject: [PATCH 21/71] better error catching --- openmc/tallies.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index b795fc83a2..faa7c3efe0 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -471,9 +471,12 @@ class Tally(IDManagerMixin): try: import vtk import vtk.util.numpy_support as nps - except ImportError: + except ModuleNotFoundError: msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ImportError(msg) + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From fb40ba894addb3ab87226e80e631aaa6288dafa6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:18:22 +0200 Subject: [PATCH 22/71] moved comment statement --- openmc/tallies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index faa7c3efe0..bec0ddc6f7 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -467,7 +467,6 @@ class Tally(IDManagerMixin): ValueError: if no MeshFilter with appropriate mesh was found (SphericalMesh not supported) """ - # check that tally has a MeshFilter try: import vtk import vtk.util.numpy_support as nps @@ -477,6 +476,7 @@ class Tally(IDManagerMixin): except ImportError as err: raise err + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): From dacd3b75b855d8f28da8f7eedea4e6508c20fd80 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:20:25 +0200 Subject: [PATCH 23/71] moved import to voxels_to_vtk --- openmc/tallies.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index bec0ddc6f7..fbf2962b9b 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -467,15 +467,6 @@ class Tally(IDManagerMixin): ValueError: if no MeshFilter with appropriate mesh was found (SphericalMesh not supported) """ - try: - import vtk - import vtk.util.numpy_support as nps - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - # check that tally has a MeshFilter mesh = None for f in self.filters: @@ -3324,7 +3315,14 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid """ - + try: + import vtk + import vtk.util.numpy_support as nps + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err # TODO: should this be a method of Tally? vtk_grid = vtk.vtkStructuredGrid() From 9ec0a0a3df1a90783e1b0e60701f972275ffd17e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 08:42:32 +0200 Subject: [PATCH 24/71] refactoring --- openmc/tallies.py | 105 +++++++++++++++++++++++----------------------- 1 file changed, 53 insertions(+), 52 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index fbf2962b9b..a25d535e5f 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -464,61 +464,29 @@ class Tally(IDManagerMixin): filename (str): the filename (must end with .vtk) Raises: - ValueError: if no MeshFilter with appropriate mesh was found - (SphericalMesh not supported) + ValueError: if no MeshFilter was found """ + try: + import vtk + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + # check that tally has a MeshFilter mesh = None for f in self.filters: if isinstance(f, openmc.MeshFilter): - if isinstance(f.mesh, (openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh)): - mesh = f.mesh - break + mesh = f.mesh + break if not mesh: raise ValueError( - "write_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + "write_to_vtk requires a MeshFilter in the tally filters" ) - if isinstance(mesh, openmc.RegularMesh): - vtk_grid = voxels_to_vtk( - x_vals=np.linspace( - mesh.lower_left[0], - mesh.upper_right[0], - num=mesh.dimension[0] + 1, - ), - y_vals=np.linspace( - mesh.lower_left[1], - mesh.upper_right[1], - num=mesh.dimension[1] + 1, - ), - z_vals=np.linspace( - mesh.lower_left[2], - mesh.upper_right[2], - num=mesh.dimension[2] + 1, - ), - mean=self.mean, - std_dev=self.std_dev, - cylindrical=False, - ) - if isinstance(mesh, openmc.RectilinearMesh): - vtk_grid = voxels_to_vtk( - x_vals=mesh.x_grid, - y_vals=mesh.y_grid, - z_vals=mesh.z_grid, - mean=self.mean, - std_dev=self.std_dev, - cylindrical=False, - ) - elif isinstance(mesh, openmc.CylindricalMesh): - vtk_grid = voxels_to_vtk( - x_vals=mesh.r_grid, - y_vals=mesh.phi_grid, - z_vals=mesh.z_grid, - mean=self.mean, - std_dev=self.std_dev, - cylindrical=True, - ) + vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) # write the .vtk file writer = vtk.vtkStructuredGridWriter() @@ -3300,20 +3268,20 @@ class Tallies(cv.CheckedList): return cls(tallies) -def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): +def voxels_to_vtk(mesh, mean, std_dev): """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. Args: - x_vals (list): X values. - y_vals (list): Y values. - z_vals (list): Z values. + mesh (openmc.StructuredMesh): The tallied mesh (SphericalMesh is not yet supported). mean (np.array): the tally mean. std_dev (np.array): the tally standard deviation. - cylindrical (bool, optional): If set to True, cylindrical coordinates - (r, phi, z) are assumed. Defaults to True. Returns: vtkStructuredGrid: a vtk object containing tally data on the appropriate grid + + Raises: + ValueError: if mesh is not openmc.RegularMesh, openmc.RectilinearMesh or + openmc.CylindricalMesh """ try: import vtk @@ -3325,13 +3293,46 @@ def voxels_to_vtk(x_vals, y_vals, z_vals, mean, std_dev, cylindrical=True): raise err # TODO: should this be a method of Tally? + system_of_coordinates = "cartesian" + + if isinstance(mesh, openmc.RegularMesh): + print('coucou') + x_vals = np.linspace( + mesh.lower_left[0], + mesh.upper_right[0], + num=mesh.dimension[0] + 1, + ) + y_vals = np.linspace( + mesh.lower_left[1], + mesh.upper_right[1], + num=mesh.dimension[1] + 1, + ) + z_vals = np.linspace( + mesh.lower_left[2], + mesh.upper_right[2], + num=mesh.dimension[2] + 1, + ) + elif isinstance(mesh, openmc.RectilinearMesh): + x_vals = mesh.x_grid + y_vals = mesh.y_grid + z_vals = mesh.z_grid + elif isinstance(mesh, openmc.CylindricalMesh): + x_vals = mesh.r_grid + y_vals = mesh.phi_grid + z_vals = mesh.z_grid + system_of_coordinates = "cylindrical" + else: + print(type(mesh)) + raise ValueError( + "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + ) vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) # create points points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if cylindrical: # transform points to cartesian coordinates + if system_of_coordinates == "cylindrical": # transform points to cartesian coordinates points_cartesian = np.copy(points) r, phi, z = points[:, 0], points[:, 1], points[:, 2] points_cartesian[:, 0] = r * np.cos(phi) From 11cee75e46ebc7e9756d81dbfae3452ff824cc10 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:12:28 +0200 Subject: [PATCH 25/71] put write_to_vtk method at the end --- openmc/tallies.py | 74 +++++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index a25d535e5f..198f5581b2 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -457,43 +457,6 @@ class Tally(IDManagerMixin): self._std_dev = np.reshape(self._std_dev.toarray(), self.shape) self._sparse = False - def write_to_vtk(self, filename): - """Writes the tally to a vtk file - - Args: - filename (str): the filename (must end with .vtk) - - Raises: - ValueError: if no MeshFilter was found - """ - try: - import vtk - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - - # check that tally has a MeshFilter - mesh = None - for f in self.filters: - if isinstance(f, openmc.MeshFilter): - mesh = f.mesh - break - - if not mesh: - raise ValueError( - "write_to_vtk requires a MeshFilter in the tally filters" - ) - - vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - def remove_score(self, score): """Remove a score from the tally @@ -3050,6 +3013,43 @@ class Tally(IDManagerMixin): new_tally.sparse = self.sparse return new_tally + def write_to_vtk(self, filename): + """Writes the tally to a vtk file + + Args: + filename (str): the filename (must end with .vtk) + + Raises: + ValueError: if no MeshFilter was found + """ + try: + import vtk + except ModuleNotFoundError: + msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' + raise ModuleNotFoundError(msg) + except ImportError as err: + raise err + + # check that tally has a MeshFilter + mesh = None + for f in self.filters: + if isinstance(f, openmc.MeshFilter): + mesh = f.mesh + break + + if not mesh: + raise ValueError( + "write_to_vtk requires a MeshFilter in the tally filters" + ) + + vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + class Tallies(cv.CheckedList): """Collection of Tallies used for an OpenMC simulation. From bff8bc4b459030448749e587fb0fde2cf53d48f9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:16:40 +0200 Subject: [PATCH 26/71] added tests --- tests/unit_tests/test_tallies.py | 62 ++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index aeeb0612d5..b6fc6cc478 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -1,4 +1,8 @@ import numpy as np +import pytest +import vtk +from os.path import exists + import openmc @@ -38,3 +42,61 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores + + +cylinder_mesh = openmc.CylindricalMesh() +cylinder_mesh.r_grid = np.linspace(1, 2, num=30) +cylinder_mesh.phi_grid = np.linspace(0, np.pi / 2, num=50) +cylinder_mesh.z_grid = np.linspace(0, 1, num=30) + +regular_mesh = openmc.RegularMesh() +regular_mesh.lower_left = [0, 0, 0] +regular_mesh.upper_right = [1, 1, 1] +regular_mesh.dimension = [10, 5, 6] + +rectilinear_mesh = openmc.RectilinearMesh() +rectilinear_mesh.x_grid = np.linspace(0, 1) +rectilinear_mesh.y_grid = np.linspace(0, 1) +rectilinear_mesh.z_grid = np.linspace(0, 1) + + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +def test_voxels_to_vtk(mesh): + vtk_grid = openmc.voxels_to_vtk(mesh, mean=None, std_dev=None) + assert isinstance(vtk_grid, vtk.vtkStructuredGrid) + + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +def test_write_to_vtk(mesh, tmpdir): + # build + tally = openmc.Tally() + tally.filters = [openmc.MeshFilter(mesh)] + filename = tmpdir / "out.vtk" + # run + tally.write_to_vtk(filename) + # test + assert exists(filename) + + +def test_write_to_vtk_raises_error_when_no_meshfilter(): + # build + tally = openmc.Tally() + + # test + expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" + with pytest.raises(ValueError, match=expected_err_msg): + tally.write_to_vtk("out.vtk") + + +def test_voxel_to_vtk_raises_error_with_wrong_mesh(): + # build + tally = openmc.Tally() + spherical_mesh = openmc.SphericalMesh() + spherical_mesh.r_grid = np.linspace(1, 2) + spherical_mesh.phi_grid = np.linspace(1, 2) + spherical_mesh.theta_grid = np.linspace(1, 2) + tally.filters = [openmc.MeshFilter(spherical_mesh)] + # test + expected_err_msg = "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" + with pytest.raises(ValueError, match=expected_err_msg): + tally.write_to_vtk("out.vtk") \ No newline at end of file From b61ec91d781caf8e7ca593a98c76d916a8efa94e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:50:28 +0200 Subject: [PATCH 27/71] removed print statement --- openmc/tallies.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 198f5581b2..88f747249d 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3296,7 +3296,6 @@ def voxels_to_vtk(mesh, mean, std_dev): system_of_coordinates = "cartesian" if isinstance(mesh, openmc.RegularMesh): - print('coucou') x_vals = np.linspace( mesh.lower_left[0], mesh.upper_right[0], From cac84d6e1f76b69d28da8674fcf81011170cf446 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 09:50:42 +0200 Subject: [PATCH 28/71] removd print statement --- openmc/tallies.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 88f747249d..020d4a3d9e 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3321,7 +3321,6 @@ def voxels_to_vtk(mesh, mean, std_dev): z_vals = mesh.z_grid system_of_coordinates = "cylindrical" else: - print(type(mesh)) raise ValueError( "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" ) From 7be28ef156a96fb4c525130d9c52e6966cf43539 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:08:58 +0200 Subject: [PATCH 29/71] added vtk_grid to meshes --- openmc/mesh.py | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 7f7b078795..dbdfb10b02 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -621,7 +621,45 @@ class RegularMesh(StructuredMesh): root_cell.fill = lattice return root_cell, cells + + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + x_vals = np.linspace( + self.lower_left[0], + self.upper_right[0], + num=self.dimension[0] + 1, + ) + y_vals = np.linspace( + self.lower_left[1], + self.upper_right[1], + num=self.dimension[1] + 1, + ) + z_vals = np.linspace( + self.lower_left[2], + self.upper_right[2], + num=self.dimension[2] + 1, + ) + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # create points + pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid def Mesh(*args, **kwargs): warnings.warn("Mesh has been renamed RegularMesh. Future versions of " @@ -821,6 +859,33 @@ class RectilinearMesh(StructuredMesh): return element + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + + x_vals = self.x_grid + y_vals = self.y_grid + z_vals = self.z_grid + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) + + # create points + pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid + class CylindricalMesh(StructuredMesh): """A 3D cylindrical mesh @@ -1012,6 +1077,37 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) + def vtk_grid(self, filename=None): + import vtk + from vtk.util import numpy_support as nps + + r_vals = self.r_grid + phi_vals = self.phi_grid + z_vals = self.z_grid + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(r_vals), len(phi_vals), len(z_vals)) + + # create points + pts_cylindrical = np.array([[r, phi, z] for z in z_vals for phi in phi_vals for r in r_vals]) + pts_cartesian = np.copy(pts_cylindrical) + r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + pts_cartesian[:, 0] = r * np.cos(phi) + pts_cartesian[:, 1] = r * np.sin(phi) + pts_cartesian[:, 2] = z + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class SphericalMesh(StructuredMesh): """A 3D spherical mesh @@ -1204,6 +1300,9 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) + def vtk_grid(self, filename=None): + raise NotImplementedError("vtk_grid not implemented for SphericalMesh") + class UnstructuredMesh(MeshBase): """A 3D unstructured mesh From 4ab6617bc4fe681bff68bce6dd4b22bf57a821d9 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:09:18 +0200 Subject: [PATCH 30/71] write_to_vtk makes use of mesh.vtk_grid() --- openmc/tallies.py | 113 +++++----------------------------------------- 1 file changed, 12 insertions(+), 101 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 020d4a3d9e..e67365df66 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3042,7 +3042,18 @@ class Tally(IDManagerMixin): "write_to_vtk requires a MeshFilter in the tally filters" ) - vtk_grid = voxels_to_vtk(mesh, self.mean, self.std_dev) + vtk_grid = mesh.vtk_grid() + + # add mean and std dev data + mean_array = vtk.vtkDoubleArray() + mean_array.SetName("mean") + mean_array.SetArray(self.mean, self.mean.size, True) + vtk_grid.GetCellData().AddArray(mean_array) + + std_dev_array = vtk.vtkDoubleArray() + std_dev_array.SetName("std_dev") + std_dev_array.SetArray(self.std_dev, self.std_dev.size, True) + vtk_grid.GetCellData().AddArray(std_dev_array) # write the .vtk file writer = vtk.vtkStructuredGridWriter() @@ -3266,103 +3277,3 @@ class Tallies(cv.CheckedList): tallies.append(tally) return cls(tallies) - - -def voxels_to_vtk(mesh, mean, std_dev): - """Creates a vtk object from a list of X, Y, Z values and mean/std_dev data. - - Args: - mesh (openmc.StructuredMesh): The tallied mesh (SphericalMesh is not yet supported). - mean (np.array): the tally mean. - std_dev (np.array): the tally standard deviation. - - Returns: - vtkStructuredGrid: a vtk object containing tally data on the appropriate grid - - Raises: - ValueError: if mesh is not openmc.RegularMesh, openmc.RectilinearMesh or - openmc.CylindricalMesh - """ - try: - import vtk - import vtk.util.numpy_support as nps - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - # TODO: should this be a method of Tally? - - system_of_coordinates = "cartesian" - - if isinstance(mesh, openmc.RegularMesh): - x_vals = np.linspace( - mesh.lower_left[0], - mesh.upper_right[0], - num=mesh.dimension[0] + 1, - ) - y_vals = np.linspace( - mesh.lower_left[1], - mesh.upper_right[1], - num=mesh.dimension[1] + 1, - ) - z_vals = np.linspace( - mesh.lower_left[2], - mesh.upper_right[2], - num=mesh.dimension[2] + 1, - ) - elif isinstance(mesh, openmc.RectilinearMesh): - x_vals = mesh.x_grid - y_vals = mesh.y_grid - z_vals = mesh.z_grid - elif isinstance(mesh, openmc.CylindricalMesh): - x_vals = mesh.r_grid - y_vals = mesh.phi_grid - z_vals = mesh.z_grid - system_of_coordinates = "cylindrical" - else: - raise ValueError( - "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" - ) - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - - # create points - points = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - if system_of_coordinates == "cylindrical": # transform points to cartesian coordinates - points_cartesian = np.copy(points) - r, phi, z = points[:, 0], points[:, 1], points[:, 2] - points_cartesian[:, 0] = r * np.cos(phi) - points_cartesian[:, 1] = r * np.sin(phi) - points_cartesian[:, 2] = z - points = points_cartesian - - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # add mean and std dev data - mean_array = vtk.vtkDoubleArray() - mean_array.SetName("mean") - if mean is None: - mean = np.zeros( - (len(x_vals) - 1) - * (len(y_vals) - 1) - * (len(z_vals) - 1) - ) - mean_array.SetArray(mean, mean.size, True) - vtk_grid.GetCellData().AddArray(mean_array) - - std_dev_array = vtk.vtkDoubleArray() - std_dev_array.SetName("std_dev") - if std_dev is None: - std_dev = np.zeros( - (len(x_vals) - 1) - * (len(y_vals) - 1) - * (len(z_vals) - 1) - ) - std_dev_array.SetArray(std_dev, std_dev.size, True) - vtk_grid.GetCellData().AddArray(std_dev_array) - - return vtk_grid From 8958b74da73e47757147ab80b22b325926291142 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:09:32 +0200 Subject: [PATCH 31/71] adapted tests --- tests/unit_tests/test_tallies.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index b6fc6cc478..e49596f6da 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -43,6 +43,23 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores +def run_dummy_sim(tally): + mat = openmc.Material() + mat.add_nuclide('Zr90', 1.0) + mat.set_density('g/cm3', 1.0) + + model = openmc.Model() + sph = openmc.Sphere(r=25.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model.geometry = openmc.Geometry([cell]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.particles = 50 + + model.tallies = openmc.Tallies([tally]) + + model.run() cylinder_mesh = openmc.CylindricalMesh() cylinder_mesh.r_grid = np.linspace(1, 2, num=30) @@ -59,19 +76,13 @@ rectilinear_mesh.x_grid = np.linspace(0, 1) rectilinear_mesh.y_grid = np.linspace(0, 1) rectilinear_mesh.z_grid = np.linspace(0, 1) - -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) -def test_voxels_to_vtk(mesh): - vtk_grid = openmc.voxels_to_vtk(mesh, mean=None, std_dev=None) - assert isinstance(vtk_grid, vtk.vtkStructuredGrid) - - @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) def test_write_to_vtk(mesh, tmpdir): # build tally = openmc.Tally() tally.filters = [openmc.MeshFilter(mesh)] filename = tmpdir / "out.vtk" + run_dummy_sim(tally) # run tally.write_to_vtk(filename) # test @@ -97,6 +108,6 @@ def test_voxel_to_vtk_raises_error_with_wrong_mesh(): spherical_mesh.theta_grid = np.linspace(1, 2) tally.filters = [openmc.MeshFilter(spherical_mesh)] # test - expected_err_msg = "voxels_to_vtk only works with openmc.RegularMesh, openmc.RectilinearMesh, openmc.CylindricalMesh" - with pytest.raises(ValueError, match=expected_err_msg): + expected_err_msg = "vtk_grid not implemented for SphericalMesh" + with pytest.raises(NotImplementedError, match=expected_err_msg): tally.write_to_vtk("out.vtk") \ No newline at end of file From 2a78c06cf512e579eaeef007a63a2c90da9cb4ed Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:21:45 +0200 Subject: [PATCH 32/71] added docstrings --- openmc/mesh.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index dbdfb10b02..f6e3953282 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -623,6 +623,15 @@ class RegularMesh(StructuredMesh): return root_cell, cells def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -860,6 +869,15 @@ class RectilinearMesh(StructuredMesh): return element def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -1078,6 +1096,15 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ import vtk from vtk.util import numpy_support as nps @@ -1301,6 +1328,16 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) def vtk_grid(self, filename=None): + """Creates a VTK object of the mesh + + Args: + filename (str, optional): Name of the vtk file to write = + (must end with .vtk). Defaults to None. + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ + # FIXME raise NotImplementedError("vtk_grid not implemented for SphericalMesh") From 913bc99564b97490d90152aa68ed9a65838bdc3e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:33:17 +0200 Subject: [PATCH 33/71] minor refactore --- openmc/mesh.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index f6e3953282..d89a980ac6 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1108,15 +1108,12 @@ class CylindricalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps - r_vals = self.r_grid - phi_vals = self.phi_grid - z_vals = self.z_grid vtk_grid = vtk.vtkStructuredGrid() - vtk_grid.SetDimensions(len(r_vals), len(phi_vals), len(z_vals)) + vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) # create points - pts_cylindrical = np.array([[r, phi, z] for z in z_vals for phi in phi_vals for r in r_vals]) + pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] pts_cartesian[:, 0] = r * np.cos(phi) From aeb66509098b597a5cd53fca88d5e21096c691a0 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:54:52 +0200 Subject: [PATCH 34/71] implemented spherical grid --- openmc/mesh.py | 29 +++++++++++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index d89a980ac6..260ed5f43c 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1334,8 +1334,33 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object """ - # FIXME - raise NotImplementedError("vtk_grid not implemented for SphericalMesh") + import vtk + from vtk.util import numpy_support as nps + + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) + + # create points + pts_cylindrical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_cartesian = np.copy(pts_cylindrical) + r, theta, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 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) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) + vtk_grid.SetPoints(vtkPts) + + if filename: + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class UnstructuredMesh(MeshBase): From 7de718a85c1ceaa0a6c485184829357d0f5c6372 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 11:56:32 +0200 Subject: [PATCH 35/71] adapted tests --- tests/unit_tests/test_tallies.py | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index e49596f6da..e3f9a486da 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -76,7 +76,12 @@ rectilinear_mesh.x_grid = np.linspace(0, 1) rectilinear_mesh.y_grid = np.linspace(0, 1) rectilinear_mesh.z_grid = np.linspace(0, 1) -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh]) +spherical_mesh = openmc.SphericalMesh() +spherical_mesh.r_grid = np.linspace(1, 2) +spherical_mesh.phi_grid = np.linspace(1, 2) +spherical_mesh.theta_grid = np.linspace(1, 2) + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_to_vtk(mesh, tmpdir): # build tally = openmc.Tally() @@ -97,17 +102,3 @@ def test_write_to_vtk_raises_error_when_no_meshfilter(): expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" with pytest.raises(ValueError, match=expected_err_msg): tally.write_to_vtk("out.vtk") - - -def test_voxel_to_vtk_raises_error_with_wrong_mesh(): - # build - tally = openmc.Tally() - spherical_mesh = openmc.SphericalMesh() - spherical_mesh.r_grid = np.linspace(1, 2) - spherical_mesh.phi_grid = np.linspace(1, 2) - spherical_mesh.theta_grid = np.linspace(1, 2) - tally.filters = [openmc.MeshFilter(spherical_mesh)] - # test - expected_err_msg = "vtk_grid not implemented for SphericalMesh" - with pytest.raises(NotImplementedError, match=expected_err_msg): - tally.write_to_vtk("out.vtk") \ No newline at end of file From 0f23bf2109505f1ebb2de692552f096b25a5e93b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:21:59 +0200 Subject: [PATCH 36/71] write_data_to_vtk SphericalMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 260ed5f43c..c780827ea6 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1324,7 +1324,7 @@ class SphericalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_t), V_p) - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -1337,6 +1337,19 @@ class SphericalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) @@ -1353,12 +1366,26 @@ class SphericalMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 15cc2a0804d3fa146bf6e337f4987ae634490513 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:23:55 +0200 Subject: [PATCH 37/71] fixed variable name --- openmc/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index c780827ea6..3dbc170d60 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1355,9 +1355,9 @@ class SphericalMesh(StructuredMesh): vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) # create points - pts_cylindrical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) - pts_cartesian = np.copy(pts_cylindrical) - r, theta, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_cartesian = np.copy(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) From 3b1d1756e4716e5e70f58a3142dc7efaed2d7913 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:26:29 +0200 Subject: [PATCH 38/71] write_data_to_vtk CylindricalMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3dbc170d60..029d714f0b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1095,7 +1095,7 @@ class CylindricalMesh(StructuredMesh): return np.multiply.outer(np.outer(V_r, V_p), V_z) - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -1108,6 +1108,19 @@ class CylindricalMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + vtk_grid = vtk.vtkStructuredGrid() vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) @@ -1124,12 +1137,26 @@ class CylindricalMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 722b572dc314ba3dc5802da57a43e0d124323556 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:31:35 +0200 Subject: [PATCH 39/71] write_data_to_vtk RegularMesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 029d714f0b..179e9a1713 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -622,7 +622,7 @@ class RegularMesh(StructuredMesh): return root_cell, cells - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -635,6 +635,19 @@ class RegularMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1] * self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1] * self.dimension[2] + cv.check_type('label', label, str) + x_vals = np.linspace( self.lower_left[0], self.upper_right[0], @@ -661,12 +674,26 @@ class RegularMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 5ecac76b6960f635f375fd48f3cb567733073d3f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:34:05 +0200 Subject: [PATCH 40/71] write_data_to_vtk rectilinear mesh --- openmc/mesh.py | 41 ++++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 179e9a1713..db2e1593c5 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -895,7 +895,7 @@ class RectilinearMesh(StructuredMesh): return element - def vtk_grid(self, filename=None): + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh Args: @@ -908,6 +908,19 @@ class RectilinearMesh(StructuredMesh): import vtk from vtk.util import numpy_support as nps + if self.volumes is None and volume_normalization: + raise RuntimeError("No volume data is present on this " + "unstructured mesh. Please load the " + " mesh information from a statepoint file.") + + # check that the data sets are appropriately sized + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + else: + assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + cv.check_type('label', label, str) + x_vals = self.x_grid y_vals = self.y_grid z_vals = self.z_grid @@ -922,12 +935,26 @@ class RectilinearMesh(StructuredMesh): vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) vtk_grid.SetPoints(vtkPts) - if filename: - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(filename) + writer.SetInputData(vtk_grid) + writer.Write() return vtk_grid From 3233ffc759708b6df7a12c8234404897fa329e51 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 12:34:58 +0200 Subject: [PATCH 41/71] removed Tally.write_to_vtk --- openmc/tallies.py | 48 ----------------------------------------------- 1 file changed, 48 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index e67365df66..d0355f14ee 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -3013,54 +3013,6 @@ class Tally(IDManagerMixin): new_tally.sparse = self.sparse return new_tally - def write_to_vtk(self, filename): - """Writes the tally to a vtk file - - Args: - filename (str): the filename (must end with .vtk) - - Raises: - ValueError: if no MeshFilter was found - """ - try: - import vtk - except ModuleNotFoundError: - msg = 'Python package vtk was not found, please install vtk to use Tally.write_to_vtk.' - raise ModuleNotFoundError(msg) - except ImportError as err: - raise err - - # check that tally has a MeshFilter - mesh = None - for f in self.filters: - if isinstance(f, openmc.MeshFilter): - mesh = f.mesh - break - - if not mesh: - raise ValueError( - "write_to_vtk requires a MeshFilter in the tally filters" - ) - - vtk_grid = mesh.vtk_grid() - - # add mean and std dev data - mean_array = vtk.vtkDoubleArray() - mean_array.SetName("mean") - mean_array.SetArray(self.mean, self.mean.size, True) - vtk_grid.GetCellData().AddArray(mean_array) - - std_dev_array = vtk.vtkDoubleArray() - std_dev_array.SetName("std_dev") - std_dev_array.SetArray(self.std_dev, self.std_dev.size, True) - vtk_grid.GetCellData().AddArray(std_dev_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - class Tallies(cv.CheckedList): """Collection of Tallies used for an OpenMC simulation. From c58465214b96f7245979f792f4bf98201551e560 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:01:47 +0200 Subject: [PATCH 42/71] removed tally tests --- tests/unit_tests/test_tallies.py | 63 -------------------------------- 1 file changed, 63 deletions(-) diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index e3f9a486da..bfff741a97 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -1,7 +1,4 @@ import numpy as np -import pytest -import vtk -from os.path import exists import openmc @@ -42,63 +39,3 @@ def test_xml_roundtrip(run_in_tmpdir): assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type assert new_tally.triggers[0].threshold == tally.triggers[0].threshold assert new_tally.triggers[0].scores == tally.triggers[0].scores - -def run_dummy_sim(tally): - mat = openmc.Material() - mat.add_nuclide('Zr90', 1.0) - mat.set_density('g/cm3', 1.0) - - model = openmc.Model() - sph = openmc.Sphere(r=25.0, boundary_type='vacuum') - cell = openmc.Cell(fill=mat, region=-sph) - model.geometry = openmc.Geometry([cell]) - - model.settings.run_mode = 'fixed source' - model.settings.batches = 2 - model.settings.particles = 50 - - model.tallies = openmc.Tallies([tally]) - - model.run() - -cylinder_mesh = openmc.CylindricalMesh() -cylinder_mesh.r_grid = np.linspace(1, 2, num=30) -cylinder_mesh.phi_grid = np.linspace(0, np.pi / 2, num=50) -cylinder_mesh.z_grid = np.linspace(0, 1, num=30) - -regular_mesh = openmc.RegularMesh() -regular_mesh.lower_left = [0, 0, 0] -regular_mesh.upper_right = [1, 1, 1] -regular_mesh.dimension = [10, 5, 6] - -rectilinear_mesh = openmc.RectilinearMesh() -rectilinear_mesh.x_grid = np.linspace(0, 1) -rectilinear_mesh.y_grid = np.linspace(0, 1) -rectilinear_mesh.z_grid = np.linspace(0, 1) - -spherical_mesh = openmc.SphericalMesh() -spherical_mesh.r_grid = np.linspace(1, 2) -spherical_mesh.phi_grid = np.linspace(1, 2) -spherical_mesh.theta_grid = np.linspace(1, 2) - -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) -def test_write_to_vtk(mesh, tmpdir): - # build - tally = openmc.Tally() - tally.filters = [openmc.MeshFilter(mesh)] - filename = tmpdir / "out.vtk" - run_dummy_sim(tally) - # run - tally.write_to_vtk(filename) - # test - assert exists(filename) - - -def test_write_to_vtk_raises_error_when_no_meshfilter(): - # build - tally = openmc.Tally() - - # test - expected_err_msg = "write_to_vtk requires a MeshFilter in the tally filters" - with pytest.raises(ValueError, match=expected_err_msg): - tally.write_to_vtk("out.vtk") From 15b98ea541fd27f6e02b5939bd482155cf44b56b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:02:12 +0200 Subject: [PATCH 43/71] added a test to meshes --- tests/unit_tests/test_mesh_to_vtk.py | 35 ++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 tests/unit_tests/test_mesh_to_vtk.py diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py new file mode 100644 index 0000000000..b5c8ca60a9 --- /dev/null +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -0,0 +1,35 @@ +import numpy as np +from os.path import exists +import pytest + +import openmc + + +regular_mesh = openmc.RegularMesh() +regular_mesh.lower_left = (0, 0, 0) +regular_mesh.upper_right = (1, 1, 1) +regular_mesh.dimension = [30, 20, 10] + +rectilinear_mesh = openmc.RectilinearMesh() +rectilinear_mesh.x_grid = np.linspace(1, 2, num=30) +rectilinear_mesh.y_grid = np.linspace(1, 2, num=30) +rectilinear_mesh.z_grid = np.linspace(1, 2, num=30) + +cylinder_mesh = openmc.CylindricalMesh() +cylinder_mesh.r_grid = np.linspace(1, 2, num=30) +cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50) +cylinder_mesh.z_grid = np.linspace(0, 1, num=30) + +spherical_mesh = openmc.SphericalMesh() +spherical_mesh.r_grid = np.linspace(1, 2, num=30) +spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50) +spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +def test_write_data_to_vtk(mesh, tmpdir): + filename = tmpdir / "out.vtk" + + data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) + + mesh.write_data_to_vtk(filename=filename, datasets={"label": data}) + assert exists(filename) \ No newline at end of file From 5174145c1ec9daacac8cfcd33e4b9212d9889fa1 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:21:36 +0200 Subject: [PATCH 44/71] added checks to test --- tests/unit_tests/test_mesh_to_vtk.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index b5c8ca60a9..e487b3668b 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,6 +1,8 @@ import numpy as np from os.path import exists import pytest +import vtk +from vtk.util import numpy_support as nps import openmc @@ -27,9 +29,30 @@ spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_data_to_vtk(mesh, tmpdir): + # BUILD filename = tmpdir / "out.vtk" data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) - mesh.write_data_to_vtk(filename=filename, datasets={"label": data}) - assert exists(filename) \ No newline at end of file + # RUN + mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) + + # TEST + assert exists(filename) + + # read file + reader = vtk.vtkStructuredGridReader() + reader.SetFileName(filename) + reader.Update() + + # check name of datasets + vtk_grid = reader.GetOutput() + array1 = vtk_grid.GetCellData().GetArray(0) + array2 = vtk_grid.GetCellData().GetArray(1) + + assert array1.GetName() == "label1" + assert array2.GetName() == "label2" + + # check size of datasets + assert nps.vtk_to_numpy(array1).size == data.size + assert nps.vtk_to_numpy(array2).size == data.size \ No newline at end of file From 5603e328b8650e874b2cf0235bba6fe792850129 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 13:24:44 +0200 Subject: [PATCH 45/71] docstrings --- openmc/mesh.py | 48 ++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index db2e1593c5..1372c5196e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -626,8 +626,16 @@ class RegularMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -899,8 +907,16 @@ class RectilinearMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -1153,8 +1169,16 @@ class CylindricalMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object @@ -1409,8 +1433,16 @@ class SphericalMesh(StructuredMesh): """Creates a VTK object of the mesh Args: - filename (str, optional): Name of the vtk file to write = - (must end with .vtk). Defaults to None. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when no volume data is found and volume_normalisation + is required Returns: vtk.vtkStructuredGrid: the VTK object From c77eb6097d6584820256ad497520273669299f9a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:18:31 +0200 Subject: [PATCH 46/71] removed check volume for structured meshes --- openmc/mesh.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1372c5196e..a63834941f 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -633,21 +633,12 @@ class RegularMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -914,21 +905,12 @@ class RectilinearMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -1176,21 +1158,12 @@ class CylindricalMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): @@ -1440,21 +1413,12 @@ class SphericalMesh(StructuredMesh): normalize the data by the volume of the mesh elements. Defaults to True. - Raises: - RuntimeError: when no volume data is found and volume_normalisation - is required - Returns: vtk.vtkStructuredGrid: the VTK object """ import vtk from vtk.util import numpy_support as nps - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") - # check that the data sets are appropriately sized for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): From 1a7c0791bf92e46475343323d7ac0b939bcb9a5a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:19:12 +0200 Subject: [PATCH 47/71] z is unchanged --- openmc/mesh.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index a63834941f..1ccb6e86e0 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1179,10 +1179,9 @@ class CylindricalMesh(StructuredMesh): # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) - r, phi, z = pts_cylindrical[:, 0], pts_cylindrical[:, 1], pts_cylindrical[:, 2] + r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) - pts_cartesian[:, 2] = z vtkPts = vtk.vtkPoints() vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) From 9867a4340abbbd2fa093acb1fd04ed1745ab6513 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:19:56 +0200 Subject: [PATCH 48/71] [skip ci] use .num_mesh_cells Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index e487b3668b..b57235edda 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -32,7 +32,7 @@ def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = tmpdir / "out.vtk" - data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) + data = np.random.random(mesh.num_mesh_cells) # RUN mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) From 7e074fe354fc60887c8c4d6540c28e1a6a338ebc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:20:17 +0200 Subject: [PATCH 49/71] [skip ci] str(filename) Co-authored-by: Paul Romano --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1372c5196e..e024710b48 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -699,7 +699,7 @@ class RegularMesh(StructuredMesh): # write the .vtk file writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) + writer.SetFileName(str(filename)) writer.SetInputData(vtk_grid) writer.Write() From 1155ffb739d8a80e5520c1013596864874a57cf5 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:21:15 +0200 Subject: [PATCH 50/71] added pathlib import --- tests/unit_tests/test_mesh_to_vtk.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index e487b3668b..cfb0407ca0 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,5 +1,6 @@ import numpy as np from os.path import exists +from pathlib import Path import pytest import vtk from vtk.util import numpy_support as nps From aeeb078344023ebb47ec7238b763c89a541a05fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:22:01 +0200 Subject: [PATCH 51/71] [skip ci] pathlib instead of os Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index b57235edda..3b7ccf9fa6 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -30,7 +30,7 @@ spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) @pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) def test_write_data_to_vtk(mesh, tmpdir): # BUILD - filename = tmpdir / "out.vtk" + filename = Path(tmpdir) / "out.vtk" data = np.random.random(mesh.num_mesh_cells) From f99e93b5839f4166972fa8e83d93de6eccb5b619 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 14:22:26 +0200 Subject: [PATCH 52/71] [skip ci] is_file() instead of exists() Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 3b7ccf9fa6..31fbd4c882 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -38,7 +38,7 @@ def test_write_data_to_vtk(mesh, tmpdir): mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) # TEST - assert exists(filename) + assert filename.is_file() # read file reader = vtk.vtkStructuredGridReader() From 1266e433c8bd3f525e6b1ea9bafef174238be09c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:23:06 +0200 Subject: [PATCH 53/71] removed unused import --- tests/unit_tests/test_mesh_to_vtk.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 572834f84f..5da62cdef4 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,5 +1,4 @@ import numpy as np -from os.path import exists from pathlib import Path import pytest import vtk From 58ccb82cf50a5ba76e622f96969ba4db6752a885 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:26:27 +0200 Subject: [PATCH 54/71] back to .dimension --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 5da62cdef4..451c869c2e 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -32,7 +32,7 @@ def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = Path(tmpdir) / "out.vtk" - data = np.random.random(mesh.num_mesh_cells) + data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) # RUN mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) From 890c3ad0d58fa59532dc143e56c0100e4540fd0d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:26:46 +0200 Subject: [PATCH 55/71] pytest.importorskip --- tests/unit_tests/test_mesh_to_vtk.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 451c869c2e..f651f5acb5 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,7 +1,8 @@ import numpy as np from pathlib import Path import pytest -import vtk + +vtk = pytest.importorskip("vtk") from vtk.util import numpy_support as nps import openmc From 48a0d15b1aa6bbbfa2f3b3c85252a9d05c923028 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:36:54 +0200 Subject: [PATCH 56/71] added error raise + new test --- openmc/mesh.py | 53 ++++++++++++++++++++++------ tests/unit_tests/test_mesh_to_vtk.py | 17 ++++++++- 2 files changed, 59 insertions(+), 11 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index b883a4a096..ada69887ad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -635,16 +635,22 @@ class RegularMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1] * self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1] * self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) x_vals = np.linspace( @@ -907,16 +913,22 @@ class RectilinearMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) x_vals = self.x_grid @@ -1160,16 +1172,22 @@ class CylindricalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) vtk_grid = vtk.vtkStructuredGrid() @@ -1414,16 +1432,23 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells """ import vtk from vtk.util import numpy_support as nps # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2] + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2] + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + cv.check_type('label', label, str) vtk_grid = vtk.vtkStructuredGrid() @@ -1642,6 +1667,11 @@ class UnstructuredMesh(MeshBase): volume_normalization : bool Whether or not to normalize the data by the volume of the mesh elements + + Raises + ------ + RuntimeError + when the size of a dataset doesn't match the number of cells """ import vtk @@ -1658,11 +1688,14 @@ class UnstructuredMesh(MeshBase): " mesh information from a statepoint file.") # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" for label, dataset in datasets.items(): if isinstance(dataset, np.ndarray): - assert dataset.size == self.n_elements + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) else: - assert len(dataset) == self.n_elements + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) cv.check_type('label', label, str) # create data arrays for the cells/points diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index f651f5acb5..9022feecd0 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -56,4 +56,19 @@ def test_write_data_to_vtk(mesh, tmpdir): # check size of datasets assert nps.vtk_to_numpy(array1).size == data.size - assert nps.vtk_to_numpy(array2).size == data.size \ No newline at end of file + assert nps.vtk_to_numpy(array2).size == data.size + +@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +def test_write_data_to_vtk_size_mismatch(mesh): + """Checks that an error is raised when the size of the dataset + doesn't match the mesh number of cells + + Args: + mesh (openmc.StructuredMesh): the mesh to test + """ + right_size = mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2] + data = np.random.random(right_size + 1) + + expected_error_msg = "The size of the dataset label should be equal to the number of cells" + with pytest.raises(RuntimeError, match=expected_error_msg): + mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) \ No newline at end of file From 8857b769324f299ad410926b2542c7b7a24a895f Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:50:28 +0200 Subject: [PATCH 57/71] refactored by adding StructuredMesh.write_data_to_vtk() --- openmc/mesh.py | 260 +++++++++++++------------------------------------ 1 file changed, 70 insertions(+), 190 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index ada69887ad..5083b215cf 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -193,6 +193,51 @@ class StructuredMesh(MeshBase): s1 = (slice(1, None),)*ndim + (slice(None),) return (vertices[s0] + vertices[s1]) / 2 + def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): + import vtk + from vtk.util import numpy_support as nps + + # check that the data sets are appropriately sized + errmsg = "The size of the dataset {} should be equal to the number of cells" + for label, dataset in datasets.items(): + if isinstance(dataset, np.ndarray): + if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + else: + if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: + raise RuntimeError(errmsg.format(label)) + cv.check_type('label', label, str) + + vtk_grid = vtk.vtkStructuredGrid() + + vtk_grid.SetDimensions(*self.dimension) + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(nps.numpy_to_vtk(points, deep=True)) + vtk_grid.SetPoints(vtkPts) + + # create VTK arrays for each of + # the data sets + for label, dataset in datasets.items(): + dataset = np.asarray(dataset).flatten() + + if volume_normalization: + dataset /= self.volumes.flatten() + + dataset_array = vtk.vtkDoubleArray() + dataset_array.SetName(label) + dataset_array.SetArray(nps.numpy_to_vtk(dataset), + dataset.size, + True) + vtk_grid.GetCellData().AddArray(dataset_array) + + # write the .vtk file + writer = vtk.vtkStructuredGridWriter() + writer.SetFileName(str(filename)) + writer.SetInputData(vtk_grid) + writer.Write() + + return vtk_grid class RegularMesh(StructuredMesh): """A regular Cartesian mesh in one, two, or three dimensions @@ -635,23 +680,7 @@ class RegularMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) x_vals = np.linspace( self.lower_left[0], @@ -668,39 +697,16 @@ class RegularMesh(StructuredMesh): self.upper_right[2], num=self.dimension[2] + 1, ) - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) # create points pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(str(filename)) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) def Mesh(*args, **kwargs): warnings.warn("Mesh has been renamed RegularMesh. Future versions of " @@ -913,60 +919,16 @@ class RectilinearMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) - - x_vals = self.x_grid - y_vals = self.y_grid - z_vals = self.z_grid - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(x_vals), len(y_vals), len(z_vals)) - # create points - pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) + pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid]) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class CylindricalMesh(StructuredMesh): @@ -1172,28 +1134,7 @@ class CylindricalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) - - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(self.r_grid), len(self.phi_grid), len(self.z_grid)) - # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) pts_cartesian = np.copy(pts_cylindrical) @@ -1201,32 +1142,12 @@ class CylindricalMesh(StructuredMesh): pts_cartesian[:, 0] = r * np.cos(phi) pts_cartesian[:, 1] = r * np.sin(phi) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class SphericalMesh(StructuredMesh): """A 3D spherical mesh @@ -1432,28 +1353,7 @@ class SphericalMesh(StructuredMesh): Returns: vtk.vtkStructuredGrid: the VTK object - - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells """ - import vtk - from vtk.util import numpy_support as nps - - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - - cv.check_type('label', label, str) - - vtk_grid = vtk.vtkStructuredGrid() - - vtk_grid.SetDimensions(len(self.r_grid), len(self.theta_grid), len(self.phi_grid)) # create points pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) @@ -1463,32 +1363,12 @@ class SphericalMesh(StructuredMesh): pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) pts_cartesian[:, 2] = r * np.cos(phi) - vtkPts = vtk.vtkPoints() - vtkPts.SetData(nps.numpy_to_vtk(pts_cartesian, deep=True)) - vtk_grid.SetPoints(vtkPts) - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() - - if volume_normalization: - dataset /= self.volumes.flatten() - - dataset_array = vtk.vtkDoubleArray() - dataset_array.SetName(label) - dataset_array.SetArray(nps.numpy_to_vtk(dataset), - dataset.size, - True) - vtk_grid.GetCellData().AddArray(dataset_array) - - # write the .vtk file - writer = vtk.vtkStructuredGridWriter() - writer.SetFileName(filename) - writer.SetInputData(vtk_grid) - writer.Write() - - return vtk_grid + return super().write_data_to_vtk( + points=pts_cartesian, + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization + ) class UnstructuredMesh(MeshBase): From 0c4de41f5076acbd93efaad6cdb970e26c890964 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 14:56:27 +0200 Subject: [PATCH 58/71] docstrings --- openmc/mesh.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 5083b215cf..449edb8d6a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -194,6 +194,24 @@ class StructuredMesh(MeshBase): return (vertices[s0] + vertices[s1]) / 2 def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): + """Creates a VTK object of the mesh + + Args: + points (list or np.array): List of (X,Y,Y) tuples. + filename (str): Name of the VTK file to write. + datasets (dict): Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization (bool, optional): Whether or not to + normalize the data by the volume of the mesh elements. + Defaults to True. + + Raises: + RuntimeError: when the size of a dataset doesn't match the number of cells + + Returns: + vtk.vtkStructuredGrid: the VTK object + """ + import vtk from vtk.util import numpy_support as nps From ad73f03a0d712311b0d5ef80f0456d2c35b61dbe Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 15:02:20 +0200 Subject: [PATCH 59/71] converted docstrings to numpy style --- openmc/mesh.py | 129 +++++++++++++++++++++++++++++++------------------ 1 file changed, 81 insertions(+), 48 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 449edb8d6a..aa338f97e1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -196,20 +196,29 @@ class StructuredMesh(MeshBase): def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - points (list or np.array): List of (X,Y,Y) tuples. - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + points : list or np.array + List of (X,Y,Y) tuples. + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Raises: - RuntimeError: when the size of a dataset doesn't match the number of cells + Raises + ------ + RuntimeError + When the size of a dataset doesn't match the number of cells - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ import vtk @@ -688,16 +697,22 @@ class RegularMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ x_vals = np.linspace( @@ -927,16 +942,22 @@ class RectilinearMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid]) @@ -1142,16 +1163,22 @@ class CylindricalMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) @@ -1361,16 +1388,22 @@ class SphericalMesh(StructuredMesh): def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh - Args: - filename (str): Name of the VTK file to write. - datasets (dict): Dictionary whose keys are the data labels - and values are the data sets. - volume_normalization (bool, optional): Whether or not to - normalize the data by the volume of the mesh elements. - Defaults to True. + Parameters + ---------- + filename : str + Name of the VTK file to write. + datasets : dict + Dictionary whose keys are the data labels + and values are the data sets. + volume_normalization : bool, optional + Whether or not to normalize the data by + the volume of the mesh elements. + Defaults to True. - Returns: - vtk.vtkStructuredGrid: the VTK object + Returns + ------- + vtk.vtkStructuredGrid + the VTK object """ # create points From a2299063f5b4061a2eb2cd762004439fadbab57c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 15:02:52 +0200 Subject: [PATCH 60/71] converted docstrings to numpy style --- tests/unit_tests/test_mesh_to_vtk.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 9022feecd0..9247540a09 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -63,8 +63,10 @@ def test_write_data_to_vtk_size_mismatch(mesh): """Checks that an error is raised when the size of the dataset doesn't match the mesh number of cells - Args: - mesh (openmc.StructuredMesh): the mesh to test + Parameters + ---------- + mesh : openmc.StructuredMesh + The mesh to test """ right_size = mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2] data = np.random.random(right_size + 1) From d80bf7bdb596a7169f0e28d42b363d1886c43448 Mon Sep 17 00:00:00 2001 From: Patrick Myers <90068356+myerspat@users.noreply.github.com> Date: Wed, 29 Jun 2022 09:25:59 -0500 Subject: [PATCH 61/71] Added comments explaining reasoning for 2D xtensor xs in photon.h Co-authored-by: Paul Romano --- include/openmc/photon.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 50ee228c44..09fb3ba01b 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -79,7 +79,9 @@ public: Tabulated1D coherent_anomalous_real_; Tabulated1D coherent_anomalous_imag_; - // Photoionization and atomic relaxation data + // Photoionization and atomic relaxation data. Subshell cross sections are + // stored separately to improve memory access pattern when calculating the + // total cross section vector shells_; xt::xtensor cross_sections_; From 2c5a845897aeb973448f31cea2d9256754262e3e Mon Sep 17 00:00:00 2001 From: Patrick Myers <90068356+myerspat@users.noreply.github.com> Date: Wed, 29 Jun 2022 09:26:36 -0500 Subject: [PATCH 62/71] Removed an include from photon.cpp Co-authored-by: Paul Romano --- src/photon.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/photon.cpp b/src/photon.cpp index 85bd435c99..a127ac2fee 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -17,7 +17,6 @@ #include "xtensor/xmath.hpp" #include "xtensor/xoperation.hpp" #include "xtensor/xslice.hpp" -#include "xtensor/xtensor_forward.hpp" #include "xtensor/xview.hpp" #include From c38ed4d8759704e829454abfc12dfa2189dd4bf1 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 29 Jun 2022 09:37:18 -0500 Subject: [PATCH 63/71] used placeholders namespace and added comment about variables initialized --- src/photon.cpp | 6 ++++-- src/physics.cpp | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index a127ac2fee..2f9bc5e3dd 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -46,6 +46,8 @@ vector> elements; PhotonInteraction::PhotonInteraction(hid_t group) { + using namespace xt::placeholders; + // Set index of element in global vector index_ = data::elements.size(); @@ -164,8 +166,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view( - cross_sections_, xt::range(shell.threshold, xt::placeholders::_), i); + auto cross_section = + xt::view(cross_sections_, xt::range(shell.threshold, _), i); cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { diff --git a/src/physics.cpp b/src/physics.cpp index cc712affd7..bcdcf7ecf9 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -347,7 +347,8 @@ void sample_photon_reaction(Particle& p) double prob_after = prob + micro.photoelectric; if (prob_after > cutoff) { - + // Get grid index, interpolation factor, and bounding subshell + // cross sections int i_grid = micro.index_grid; double f = micro.interp_factor; const auto& xs_lower = xt::row(element.cross_sections_, i_grid); From 7cba2ad03012658481c6edbc04f070892e90dd3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:46:10 +0200 Subject: [PATCH 64/71] [skip ci] X Y Y typo Co-authored-by: Paul Romano --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index aa338f97e1..90f4c55a9a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -199,7 +199,7 @@ class StructuredMesh(MeshBase): Parameters ---------- points : list or np.array - List of (X,Y,Y) tuples. + List of (X,Y,Z) tuples. filename : str Name of the VTK file to write. datasets : dict From 5180ab3f8a6031cfbcd708c04e7fd23f8936085a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:46:43 +0200 Subject: [PATCH 65/71] [skip ci] removed default value Co-authored-by: Paul Romano --- openmc/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 90f4c55a9a..73a49bfc19 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -208,7 +208,6 @@ class StructuredMesh(MeshBase): volume_normalization : bool, optional Whether or not to normalize the data by the volume of the mesh elements. - Defaults to True. Raises ------ From 94ef0c40abe6f4a208aaaccad36c17ff486dd9b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Delaporte-Mathurin?= <40028739+RemDelaporteMathurin@users.noreply.github.com> Date: Wed, 29 Jun 2022 18:47:24 +0200 Subject: [PATCH 66/71] [skip ci] refactore upperright and lowerleft Co-authored-by: Paul Romano --- openmc/mesh.py | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 73a49bfc19..7155b0ba04 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -714,21 +714,10 @@ class RegularMesh(StructuredMesh): the VTK object """ - x_vals = np.linspace( - self.lower_left[0], - self.upper_right[0], - num=self.dimension[0] + 1, - ) - y_vals = np.linspace( - self.lower_left[1], - self.upper_right[1], - num=self.dimension[1] + 1, - ) - z_vals = np.linspace( - self.lower_left[2], - self.upper_right[2], - num=self.dimension[2] + 1, - ) + ll, ur = self.lower_left, self.upper_right + x_vals = np.linspace(ll[0], ur[0], num=self.dimension[0] + 1) + y_vals = np.linspace(ll[1], ur[1], num=self.dimension[1] + 1) + z_vals = np.linspace(ll[2], ur[2], num=self.dimension[2] + 1) # create points pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals]) From 96c89f138daa97f78cb996449a32a88ca58754ff Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 19:01:39 +0200 Subject: [PATCH 67/71] multiple lines --- openmc/mesh.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 7155b0ba04..209002e42a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1169,7 +1169,14 @@ class CylindricalMesh(StructuredMesh): the VTK object """ # create points - pts_cylindrical = np.array([[r, phi, z] for z in self.z_grid for phi in self.phi_grid for r in self.r_grid]) + pts_cylindrical = np.array( + [ + [r, phi, z] + for z in self.z_grid + for phi in self.phi_grid + for r in self.r_grid + ] + ) pts_cartesian = np.copy(pts_cylindrical) r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1] pts_cartesian[:, 0] = r * np.cos(phi) @@ -1395,7 +1402,14 @@ class SphericalMesh(StructuredMesh): """ # create points - pts_spherical = np.array([[r, theta, phi] for phi in self.phi_grid for theta in self.theta_grid for r in self.r_grid]) + pts_spherical = np.array( + [ + [r, theta, phi] + for phi in self.phi_grid + for theta in self.theta_grid + for r in self.r_grid + ] + ) pts_cartesian = np.copy(pts_spherical) r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) From fabe3e607c92891ec740aefa632eba010428eef2 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Jun 2022 19:02:25 +0200 Subject: [PATCH 68/71] pathlib --- openmc/mesh.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 209002e42a..2a53c8a91e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -698,7 +698,7 @@ class RegularMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -932,7 +932,7 @@ class RectilinearMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1153,7 +1153,7 @@ class CylindricalMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1385,7 +1385,7 @@ class SphericalMesh(StructuredMesh): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels @@ -1592,7 +1592,7 @@ class UnstructuredMesh(MeshBase): Parameters ---------- - filename : str + filename : str or pathlib.Path Name of the VTK file to write. datasets : dict Dictionary whose keys are the data labels From 0b0265cb592685d8e6daec688229996d7b7e9f5a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 1 Jul 2022 11:37:16 +0200 Subject: [PATCH 69/71] adapted tests --- tests/unit_tests/test_mesh_to_vtk.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 9247540a09..26a393b411 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -33,7 +33,7 @@ def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = Path(tmpdir) / "out.vtk" - data = np.random.random(mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2]) + data = np.random.random(mesh.num_mesh_cells) # RUN mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data}) @@ -68,7 +68,7 @@ def test_write_data_to_vtk_size_mismatch(mesh): mesh : openmc.StructuredMesh The mesh to test """ - right_size = mesh.dimension[0]*mesh.dimension[1]*mesh.dimension[2] + right_size = mesh.num_mesh_cells data = np.random.random(right_size + 1) expected_error_msg = "The size of the dataset label should be equal to the number of cells" From 9c892f04dbd36bc423a0560b72bb2da3a1b2a21a Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 1 Jul 2022 09:42:35 +0000 Subject: [PATCH 70/71] added num_mesh_cells property to StructuredMesh --- openmc/mesh.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2a53c8a91e..925a63aeec 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -193,6 +193,10 @@ class StructuredMesh(MeshBase): s1 = (slice(1, None),)*ndim + (slice(None),) return (vertices[s0] + vertices[s1]) / 2 + @property + def num_mesh_cells(self): + return np.prod(self.dimension) + def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh @@ -1547,6 +1551,7 @@ class UnstructuredMesh(MeshBase): "been loaded from a statepoint file.") return len(self._centroids) + @centroids.setter def centroids(self, centroids): cv.check_type("Unstructured mesh centroids", centroids, From 52d59a4c1a924039ca6eb790a1786fe85fce7574 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Fri, 1 Jul 2022 09:43:44 +0000 Subject: [PATCH 71/71] minor refactoring --- openmc/mesh.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 925a63aeec..d5a460baa8 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -348,10 +348,6 @@ class RegularMesh(StructuredMesh): dims = self._dimension return [(u - l) / d for u, l, d in zip(us, ls, dims)] - @property - def num_mesh_cells(self): - return np.prod(self._dimension) - @property def volumes(self): """Return Volumes for every mesh cell