diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 09b1b11442..89d973a2ba 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -139,6 +139,8 @@ The current version of the statepoint file format is 17.0. - **internal** (*int*) -- Flag indicating the presence of tally data (0) or absence of tally data (1). All user defined tallies will have a value of 0 unless otherwise instructed. + - **multiply_density** (*int*) -- Flag indicating whether reaction + rates should be multiplied by atom density (1) or not (0). :Datasets: - **n_realizations** (*int*) -- Number of realizations. - **n_filters** (*int*) -- Number of filters used. diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 897941ad40..b88877b838 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -69,6 +69,12 @@ The ```` element accepts the following sub-elements: list of valid scores can be found in the :ref:`user's guide `. + :multiply_density: + A boolean that indicates whether reaction rate scores should be computed by + multiplying by the atom density of a nuclide present in a material. + + *Default*: true + :trigger: Precision trigger applied to all filter bins and nuclides for this tally. It must specify the trigger's type, threshold and scores to which it will diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index b05c76bbb8..1cc9d29728 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -29,6 +29,7 @@ namespace openmc { class Nuclide { public: + //============================================================================ // Types, aliases using EmissionMode = ReactionProduct::EmissionMode; struct EnergyGrid { @@ -36,18 +37,32 @@ public: vector energy; }; + //============================================================================ // Constructors/destructors Nuclide(hid_t group, const vector& temperature); ~Nuclide(); + //============================================================================ + // Methods + //! Initialize logarithmic grid for energy searches void init_grid(); + //! Calculate microscopic cross sections + // + //! \param[in] i_sab Index in data::thermal_scatt + //! \param[in] i_log_union Log-grid search index + //! \param[in] sab_frac S(a,b) table fraction + //! \param[in,out] p Particle object void calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p); + //! Calculate thermal scattering cross section + // + //! \param[in] i_sab Index in data::thermal_scatt + //! \param[in] sab_frac S(a,b) table fraction + //! \param[in,out] p Particle object void calculate_sab_xs(int i_sab, double sab_frac, Particle& p); - // Methods double nu(double E, EmissionMode mode, int group = 0) const; void calculate_elastic_xs(Particle& p) const; @@ -69,6 +84,7 @@ public: double collapse_rate(int MT, double temperature, gsl::span energy, gsl::span flux) const; + //============================================================================ // Data members std::string name_; //!< Name of nuclide, e.g. "U235" int Z_; //!< Atomic number diff --git a/include/openmc/particle.h b/include/openmc/particle.h index da2f6a61fc..e3f23fd5b1 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -35,6 +35,9 @@ public: Particle() = default; + //========================================================================== + // Methods + double speed() const; //! create a secondary particle @@ -106,6 +109,16 @@ public: //! create a particle restart HDF5 file void write_restart() const; + + //! Update microscopic cross section cache + // + //! \param[in] i_nuclide Index in data::nuclides + //! \param[in] i_grid Index on log union grid + //! \param[in] i_sab Index in data::thermal_scatt + //! \param[in] sab_frac S(a,b) table fraction + //! \param[in] ncrystal_xs Thermal scattering xs from NCrystal + void update_neutron_xs(int i_nuclide, int i_grid = C_NONE, int i_sab = C_NONE, + double sab_frac = 0.0, double ncrystal_xs = -1.0); }; //============================================================================ diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index d3d00571a8..c4f959666f 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -218,8 +218,9 @@ struct BoundaryInfo { * https://doi.org/10.1016/j.anucene.2017.11.032. */ class ParticleData { - public: + //---------------------------------------------------------------------------- + // Constructors ParticleData(); private: diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index e840a0aa14..6ad3bef6f3 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -37,6 +37,8 @@ public: void set_active(bool active) { active_ = active; } + void set_multiply_density(bool value) { multiply_density_ = value; } + void set_writable(bool writable) { writable_ = writable; } void set_scores(pugi::xml_node node); @@ -62,6 +64,8 @@ public: int32_t n_filter_bins() const { return n_filter_bins_; } + bool multiply_density() const { return multiply_density_; } + bool writable() const { return writable_; } //---------------------------------------------------------------------------- @@ -139,6 +143,9 @@ private: int32_t n_filter_bins_ {0}; + //! Whether to multiply by atom density for reaction rates + bool multiply_density_ {true}; + gsl::index index_; }; diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index 85af8ffff4..83c28528e4 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -39,6 +39,9 @@ _dll.openmc_tally_get_filters.argtypes = [ c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)] _dll.openmc_tally_get_filters.restype = c_int _dll.openmc_tally_get_filters.errcheck = _error_handler +_dll.openmc_tally_get_multiply_density.argtypes = [c_int32, POINTER(c_bool)] +_dll.openmc_tally_get_multiply_density.restype = c_int +_dll.openmc_tally_get_multiply_density.errcheck = _error_handler _dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_tally_get_n_realizations.restype = c_int _dll.openmc_tally_get_n_realizations.errcheck = _error_handler @@ -75,6 +78,9 @@ _dll.openmc_tally_set_estimator.errcheck = _error_handler _dll.openmc_tally_set_id.argtypes = [c_int32, c_int32] _dll.openmc_tally_set_id.restype = c_int _dll.openmc_tally_set_id.errcheck = _error_handler +_dll.openmc_tally_set_multiply_density.argtypes = [c_int32, c_bool] +_dll.openmc_tally_set_multiply_density.restype = c_int +_dll.openmc_tally_set_multiply_density.errcheck = _error_handler _dll.openmc_tally_set_nuclides.argtypes = [c_int32, c_int, POINTER(c_char_p)] _dll.openmc_tally_set_nuclides.restype = c_int _dll.openmc_tally_set_nuclides.errcheck = _error_handler @@ -174,6 +180,10 @@ class Tally(_FortranObjectWithID): List of tally filters mean : numpy.ndarray An array containing the sample mean for each bin + multiply_density : bool + Whether reaction rates should be multiplied by atom density + + .. versionadded:: 0.13.4 nuclides : list of str List of nuclides to score results for num_realizations : int @@ -363,6 +373,16 @@ class Tally(_FortranObjectWithID): def writable(self, writable): _dll.openmc_tally_set_writable(self._index, writable) + @property + def multiply_density(self): + multiply_density = c_bool() + _dll.openmc_tally_get_multiply_density(self._index, multiply_density) + return multiply_density.value + + @multiply_density.setter + def multiply_density(self, multiply_density): + _dll.openmc_tally_set_multiply_density(self._index, multiply_density) + def reset(self): """Reset results and num_realizations of tally""" _dll.openmc_tally_reset(self._index) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 138c0c718d..7c1ecd1377 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -398,6 +398,10 @@ class StatePoint: tally._sp_filename = self._f.filename tally.name = group['name'][()].decode() if 'name' in group else '' + # Check if tally has multiply_density attribute + if "multiply_density" in group.attrs: + tally.multiply_density = group.attrs["multiply_density"].item() > 0 + # Read the number of realizations n_realizations = group['n_realizations'][()] diff --git a/openmc/tallies.py b/openmc/tallies.py index 596b324bb3..85daaa17a1 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -14,7 +14,7 @@ import scipy.sparse as sps import openmc import openmc.checkvalue as cv -from ._xml import clean_indentation, reorder_attributes +from ._xml import clean_indentation, reorder_attributes, get_text from .mixin import IDManagerMixin from .mesh import MeshBase @@ -54,6 +54,10 @@ class Tally(IDManagerMixin): Unique identifier for the tally name : str Name of the tally + multiply_density : bool + Whether reaction rates should be multiplied by atom density + + .. versionadded:: 0.13.4 filters : list of openmc.Filter List of specified filters for the tally nuclides : list of str @@ -111,6 +115,7 @@ class Tally(IDManagerMixin): self._estimator = None self._triggers = cv.CheckedList(openmc.Trigger, 'tally triggers') self._derivative = None + self._multiply_density = True self._num_realizations = 0 self._with_summary = False @@ -138,12 +143,17 @@ class Tally(IDManagerMixin): parts.append('{: <15}=\t{}'.format('Nuclides', nuclides)) parts.append('{: <15}=\t{}'.format('Scores', self.scores)) parts.append('{: <15}=\t{}'.format('Estimator', self.estimator)) + parts.append('{: <15}=\t{}'.format('Multiply dens.', self.multiply_density)) return '\n\t'.join(parts) @property def name(self): return self._name + @property + def multiply_density(self): + return self._multiply_density + @property def filters(self): return self._filters @@ -323,6 +333,11 @@ class Tally(IDManagerMixin): cv.check_type('tally name', name, str, none_ok=True) self._name = name + @multiply_density.setter + def multiply_density(self, value): + cv.check_type('multiply density', value, bool) + self._multiply_density = value + @derivative.setter def derivative(self, deriv): cv.check_type('tally derivative', deriv, openmc.TallyDerivative, @@ -829,6 +844,10 @@ class Tally(IDManagerMixin): if self.name != '': element.set("name", self.name) + # Multiply by density + if not self.multiply_density: + element.set("multiply_density", str(self.multiply_density).lower()) + # Optional Tally filters if len(self.filters) > 0: subelement = ET.SubElement(element, "filters") @@ -885,6 +904,10 @@ class Tally(IDManagerMixin): name = elem.get('name', '') tally = cls(tally_id=tally_id, name=name) + text = get_text(elem, 'multiply_density') + if text is not None: + tally.multiply_density = text in ('true', '1') + # Read filters filters_elem = elem.find('filters') if filters_elem is not None: diff --git a/src/material.cpp b/src/material.cpp index 82dd0d4d9e..610b54add2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -353,7 +353,7 @@ Material::~Material() model::material_map.erase(id_); } -Material & Material::clone() +Material& Material::clone() { std::unique_ptr mat = std::make_unique(); @@ -868,21 +868,12 @@ void Material::calculate_neutron_xs(Particle& p) const // ====================================================================== // CALCULATE MICROSCOPIC CROSS SECTION - // Determine microscopic cross sections for this nuclide + // Get nuclide index int i_nuclide = nuclide_[i]; - // Calculate microscopic cross section for this nuclide - auto& micro {p.neutron_xs(i_nuclide)}; - if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT || - i_sab != micro.index_sab || sab_frac != micro.sab_frac) { - data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p); - - // If NCrystal is being used, update micro cross section cache - if (ncrystal_xs >= 0.0) { - data::nuclides[i_nuclide]->calculate_elastic_xs(p); - ncrystal_update_micro(ncrystal_xs, micro); - } - } + // Update microscopic cross section for this nuclide + p.update_neutron_xs(i_nuclide, i_grid, i_sab, sab_frac, ncrystal_xs); + auto& micro = p.neutron_xs(i_nuclide); // ====================================================================== // ADD TO MACROSCOPIC CROSS SECTION diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 99ed44b112..8f73a53093 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -471,7 +471,7 @@ void Nuclide::create_derived( xs_cdf_sum += (std::sqrt(E[i]) * xs[i] + std::sqrt(E[i + 1]) * xs[i + 1]) / 2.0 * (E[i + 1] - E[i]); - xs_cdf_[i+1] = xs_cdf_sum; + xs_cdf_[i + 1] = xs_cdf_sum; } } } diff --git a/src/particle.cpp b/src/particle.cpp index 1ec95a4ae6..f455ab8349 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -36,6 +36,10 @@ namespace openmc { +//============================================================================== +// Particle implementation +//============================================================================== + double Particle::speed() const { // Determine mass in eV/c^2 @@ -433,7 +437,9 @@ void Particle::cross_surface() #ifdef DAGMC // in DAGMC, we know what the next cell should be if (surf->geom_type_ == GeometryType::DAG) { - int32_t i_cell = next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) - 1; + int32_t i_cell = + next_cell(i_surface, cell_last(n_coord() - 1), lowest_coord().universe) - + 1; // save material and temp material_last() = material(); sqrtkT_last() = sqrtkT(); @@ -703,6 +709,29 @@ void Particle::write_restart() const } // #pragma omp critical } +void Particle::update_neutron_xs( + int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs) +{ + // Get microscopic cross section cache + auto& micro = this->neutron_xs(i_nuclide); + + // If the cache doesn't match, recalculate micro xs + if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT || + i_sab != micro.index_sab || sab_frac != micro.sab_frac) { + data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this); + + // If NCrystal is being used, update micro cross section cache + if (ncrystal_xs >= 0.0) { + data::nuclides[i_nuclide]->calculate_elastic_xs(*this); + ncrystal_update_micro(ncrystal_xs, micro); + } + } +} + +//============================================================================== +// Non-method functions +//============================================================================== + std::string particle_type_to_str(ParticleType type) { switch (type) { diff --git a/src/state_point.cpp b/src/state_point.cpp index eefce8fa47..bff6213206 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -193,6 +193,12 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) continue; } + if (tally->multiply_density()) { + write_attribute(tally_group, "multiply_density", 1); + } else { + write_attribute(tally_group, "multiply_density", 0); + } + if (tally->estimator_ == TallyEstimator::ANALOG) { write_dataset(tally_group, "estimator", "analog"); } else if (tally->estimator_ == TallyEstimator::TRACKLENGTH) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index e29384772f..85d999fd4a 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -92,6 +92,10 @@ Tally::Tally(pugi::xml_node node) if (check_for_node(node, "name")) name_ = get_node_value(node, "name"); + if (check_for_node(node, "multiply_density")) { + multiply_density_ = get_node_value_bool(node, "multiply_density"); + } + // ======================================================================= // READ DATA FOR FILTERS @@ -564,11 +568,12 @@ void Tally::set_nuclides(const vector& nuclides) nuclides_.push_back(-1); } else { auto search = data::nuclide_map.find(nuc); - if (search == data::nuclide_map.end()) - fatal_error(fmt::format("Could not find the nuclide {} specified in " - "tally {} in any material", - nuc, id_)); - nuclides_.push_back(search->second); + if (search == data::nuclide_map.end()) { + int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0); + if (err < 0) + throw std::runtime_error {openmc_err_msg}; + } + nuclides_.push_back(data::nuclide_map.at(nuc)); } } } @@ -1108,6 +1113,28 @@ extern "C" int openmc_tally_set_writable(int32_t index, bool writable) return 0; } +extern "C" int openmc_tally_get_multiply_density(int32_t index, bool* value) +{ + if (index < 0 || index >= model::tallies.size()) { + set_errmsg("Index in tallies array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + *value = model::tallies[index]->multiply_density(); + + return 0; +} + +extern "C" int openmc_tally_set_multiply_density(int32_t index, bool value) +{ + if (index < 0 || index >= model::tallies.size()) { + set_errmsg("Index in tallies array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + model::tallies[index]->set_multiply_density(value); + + return 0; +} + extern "C" int openmc_tally_get_scores(int32_t index, int** scores, int* n) { if (index < 0 || index >= model::tallies.size()) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 68eb344af4..3e58fec05b 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -939,21 +939,17 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (i_nuclide >= 0) { const auto& micro = p.photon_xs(i_nuclide); - double xs = (score_bin == COHERENT) - ? micro.coherent - : (score_bin == INCOHERENT) ? micro.incoherent - : (score_bin == PHOTOELECTRIC) - ? micro.photoelectric - : micro.pair_production; + double xs = (score_bin == COHERENT) ? micro.coherent + : (score_bin == INCOHERENT) ? micro.incoherent + : (score_bin == PHOTOELECTRIC) ? micro.photoelectric + : micro.pair_production; score = xs * atom_density * flux; } else { - double xs = (score_bin == COHERENT) - ? p.macro_xs().coherent - : (score_bin == INCOHERENT) - ? p.macro_xs().incoherent - : (score_bin == PHOTOELECTRIC) - ? p.macro_xs().photoelectric - : p.macro_xs().pair_production; + double xs = (score_bin == COHERENT) ? p.macro_xs().coherent + : (score_bin == INCOHERENT) ? p.macro_xs().incoherent + : (score_bin == PHOTOELECTRIC) + ? p.macro_xs().photoelectric + : p.macro_xs().pair_production; score = xs * flux; } break; @@ -2308,6 +2304,9 @@ void score_tracklength_tally(Particle& p, double distance) // Determine the tracklength estimate of the flux double flux = p.wgt() * distance; + // Set 'none' value for log union grid index + int i_log_union = C_NONE; + for (auto i_tally : model::active_tracklength_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2331,11 +2330,25 @@ void score_tracklength_tally(Particle& p, double distance) double atom_density = 0.; if (i_nuclide >= 0) { if (p.material() != MATERIAL_VOID) { - auto j = - model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; - if (j == C_NONE) - continue; - atom_density = model::materials[p.material()]->atom_density_(j); + const auto& mat = model::materials[p.material()]; + auto j = mat->mat_nuclide_index_[i_nuclide]; + if (j == C_NONE) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = static_cast(ParticleType::neutron); + i_log_union = std::log(p.E() / data::energy_min[neutron]) / + simulation::log_spacing; + } + + // Update micro xs cache + if (!tally.multiply_density()) { + p.update_neutron_xs(i_nuclide, i_log_union); + atom_density = 1.0; + } + } else { + atom_density = + tally.multiply_density() ? mat->atom_density_(j) : 1.0; + } } } @@ -2371,6 +2384,9 @@ void score_collision_tally(Particle& p) flux = p.wgt_last() / p.macro_xs().total; } + // Set 'none value for log union grid index + int i_log_union = C_NONE; + for (auto i_tally : model::active_collision_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2393,11 +2409,25 @@ void score_collision_tally(Particle& p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = - model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; - if (j == C_NONE) - continue; - atom_density = model::materials[p.material()]->atom_density_(j); + const auto& mat = model::materials[p.material()]; + auto j = mat->mat_nuclide_index_[i_nuclide]; + if (j == C_NONE) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = static_cast(ParticleType::neutron); + i_log_union = std::log(p.E() / data::energy_min[neutron]) / + simulation::log_spacing; + } + + // Update micro xs cache + if (!tally.multiply_density()) { + p.update_neutron_xs(i_nuclide, i_log_union); + atom_density = 1.0; + } + } else { + atom_density = + tally.multiply_density() ? mat->atom_density_(j) : 1.0; + } } // TODO: consider replacing this "if" with pointers or templates diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat index 6a58c1725f..673d2143be 100644 --- a/tests/regression_tests/tallies/results_true.dat +++ b/tests/regression_tests/tallies/results_true.dat @@ -1 +1 @@ -6b1d8d6f4d7a70af6c39cc76b9267a61ba9a9d0c14b75a4df6f83e72a2bfaf1533caaf936c2185f0a158443eb9da266bd5a77b7996020fd638551ea841d821e0 \ No newline at end of file +d1decdbec6cb59df91ba5c42cb37a04f413a34fa7faf7ad1eecfd7d53a14af8cb65b58335c4ede4e88f6d9ab35a1746251983cc991d74eab3e678887c84183bd \ No newline at end of file diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 0633ceb183..b4dd9197b1 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -366,6 +366,19 @@ def test_tally_activate(lib_simulation_init): assert t.active +def test_tally_multiply_density(lib_simulation_init): + # multiply_density is True by default + t = openmc.lib.tallies[1] + assert t.multiply_density + + # Make sure setting multiply_density works + t.multiply_density = False + assert not t.multiply_density + + # Reset to True + t.multiply_density = True + + def test_tally_writable(lib_simulation_init): t = openmc.lib.tallies[1] assert t.writable diff --git a/tests/unit_tests/test_tally_multiply_density.py b/tests/unit_tests/test_tally_multiply_density.py new file mode 100644 index 0000000000..66ef41e0da --- /dev/null +++ b/tests/unit_tests/test_tally_multiply_density.py @@ -0,0 +1,50 @@ +import numpy as np +import openmc +import pytest + + +def test_micro_macro_compare(): + # Create simple sphere model with H1 and H2 + mat = openmc.Material() + mat.add_components({'H1': 1.0, 'H2': 1.0}) + mat.set_density('g/cm3', 1.0) + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 10 + + # Set up two reaction rate tallies, one that multplies by density and the + # other that doesn't + tally_macro = openmc.Tally() + tally_macro.nuclides = ['H1', 'H2', 'H3'] + tally_macro.scores = ['total', 'elastic'] + tally_micro = openmc.Tally() + tally_micro.nuclides = ['H1', 'H2', 'H3'] + tally_micro.scores = ['total', 'elastic'] + tally_micro.multiply_density = False + model.tallies = [tally_macro, tally_micro] + + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + tally_macro = sp.tallies[tally_macro.id] + tally_micro = sp.tallies[tally_micro.id] + + # Make sure multply_density attribute from statepoint is set correctly + assert tally_macro.multiply_density + assert not tally_micro.multiply_density + + # Dividing macro by density should give micro + density = mat.get_nuclide_atom_densities() + for nuc in ('H1', 'H2'): + micro_derived = tally_macro.get_values(nuclides=[nuc]) / density[nuc] + micro = tally_micro.get_values(nuclides=[nuc]) + assert micro_derived == pytest.approx(micro) + + # For macro tally, H3 scores should be zero + assert np.all(tally_macro.get_values(nuclides=['H3']) == 0.0) + + # For micro tally, H3 scores should be positive + assert np.all(tally_micro.get_values(nuclides=['H3']) > 0.0)