diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index db72984a63b..b25e7d517f5 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1166,11 +1166,12 @@ attributes/sub-elements: The ```` element triggers OpenMC to bank particles crossing certain surfaces and write out the source bank in a separate file called -``surface_source.h5``. One or multiple surface IDs and one cell ID can be used -to select the surfaces of interest. If no surface IDs are declared, every surface -of the model is eligible to bank particles. In that case, a cell ID (using -either the ``cell``, ``cellfrom`` or ``cellto`` attributes) can be used to select -every surface of a specific cell. This element has the following +``surface_source.h5``. One or multiple surface IDs and one or multiple cell IDs can be used +to select the surfaces of interest. The cell IDs can have direction flags assosciated with them +to further filter the banked particles. Allowed directions are ``to``, ``from`` or ``both``. +If only one cell ID is used in banking, the ``cell``, ``cellfrom`` or ``cellto`` attributes +can be used instead of the ``cells`` and ``directions`` attributes. If no surface IDs are declared, +every surface of the model is eligible to bank particles. This element has the following attributes/sub-elements: :surface_ids: @@ -1206,6 +1207,18 @@ attributes/sub-elements: .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf + :cells: + A list of integers representing the cell IDs used to determine if particles crossing + identified surfaces are to be banked. + + *Default*: None + + :directions: + A list of strings representing the directions corresponding to the cell IDs. Allowed values are + ``to``, ``from`` or ``both``. Must have the same length as ``cells``. + + *Default*: None + :cell: An integer representing the cell ID used to determine if particles crossing identified surfaces are to be banked. Particles coming from or going to this @@ -1227,9 +1240,15 @@ attributes/sub-elements: *Default*: None +.. note:: The ``cell``, ``cellfrom`` or ``cellto`` attributes cannot be + used simultaneously with ``cells`` attribute. + .. note:: The ``cell``, ``cellfrom`` and ``cellto`` attributes cannot be used simultaneously. +.. note:: If ``cells`` attribute is defined and ``directions`` attribute is not, + all the directions of in ``cells`` will default to ``both``. + .. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum" are not eligible to store any particles when using ``cell``, ``cellfrom`` or ``cellto`` attributes. It is recommended to use surface IDs instead. diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5a04fedd70a..309a6e1b582 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -323,8 +323,32 @@ crossing any surface of the model will be banked:: settings.surf_source_write = {'max_particles': 10000} -A cell ID can also be used to bank particles that are crossing any surface of -a cell that particles are either coming from or going to:: +A list of cell IDs can also be used to bank particles that are crossing surfaces of +cells that particles are either coming from or going to:: + + settings.surf_source_write = { + 'surfaces_ids': [1, 2, 3], + 'cells': [1, 2] + 'max_particles': 10000 + } + +In this example, particles that are crossing surfaces with IDs of 1, 2, or 3 and +entering or exiting cells with IDs 1 or 2 are banked. + +To account specifically for particles leaving or entering a given cell, +a list of directions can be used to further filter particles banked:: + + settings.surf_source_write = { + 'surfaces_ids': [1, 2, 3], + 'cells': [1, 2] + 'directions': '"to", "from"' + 'max_particles': 10000 + } + +In this example, particles that are crossing surfaces with IDs of 1, 2, or 3 and +entering cell with ID 1 or exiting cell with ID 2 are banked. + +If only one cell is used to filter particles, an alternative syntax can be used to bank particles:: settings.surf_source_write = {'cell': 1, 'max_particles': 10000} @@ -333,10 +357,10 @@ be banked excluding any surface that does not use a 'transmission' or 'vacuum' boundary condition. .. note:: Surfaces with boundary conditions that are not "transmission" or "vacuum" - are not eligible to store any particles when using ``cell``, ``cellfrom`` + are not eligible to store any particles when using ``cells``, ``cell``, ``cellfrom`` or ``cellto`` attributes. It is recommended to use surface IDs instead. -Surface IDs can be used in combination with a cell ID:: +Another example that combines surface IDs with a cell ID:: settings.surf_source_write = { 'cell': 1, diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..c89a218d636 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -6,6 +6,7 @@ #include #include +#include #include #include "pugixml.hpp" @@ -178,9 +179,11 @@ extern int64_t ssw_max_particles; //!< maximum number of particles to be extern int64_t ssw_max_files; //!< maximum number of surface source files //!< to be created extern int64_t ssw_cell_id; //!< Cell id for the surface source - //!< write setting -extern SSWCellType ssw_cell_type; //!< Type of option for the cell - //!< argument of surface source write + +extern std::unordered_map + ssw_cells; //!< Cell ids and directions + //!< for the surface source write setting + extern TemperatureMethod temperature_method; //!< method for choosing temperatures extern double @@ -213,6 +216,8 @@ void read_settings_xml(pugi::xml_node root); void free_memory_settings(); +SSWCellType ssw_cell_type_from_string(std::string_view s); + } // namespace openmc #endif // OPENMC_SETTINGS_H diff --git a/openmc/settings.py b/openmc/settings.py index 2ad1d06e757..e19d057f6be 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -280,15 +280,23 @@ class Settings: process (int) :max_source_files: Maximum number of surface source files to be created (int) :mcpl: Output in the form of an MCPL-file (bool) - :cell: Cell ID used to determine if particles crossing identified + :cells: List of cell IDs used to determine if particles crossing identified surfaces are to be banked. Particles coming from or going to this declared cell will be banked (int) + :directions: List of directions corresponding to cells. Acceptable entries are: + "from", "to", or "both" (str) + :cell: Cell ID used to determine if particles crossing identified + surfaces are to be banked. Particles coming from or going to this + declared cell will be banked (int) ("cell" will be deprecated in the future, + use "cells" instead.) :cellfrom: Cell ID used to determine if particles crossing identified surfaces are to be banked. Particles coming from this - declared cell will be banked (int) + declared cell will be banked (int) ("cellfrom" will be deprecated in the future, + use "cells" and "directions" instead.) :cellto: Cell ID used to determine if particles crossing identified surfaces are to be banked. Particles going to this declared - cell will be banked (int) + cell will be banked (int) ("cellto" will be deprecated in the future, + use "cells" and "directions" instead.) survival_biasing : bool Indicate whether survival biasing is to be used tabular_legendre : dict @@ -848,16 +856,22 @@ def surf_source_write(self, surf_source_write: dict): "surface source writing key", key, ("surface_ids", "max_particles", "max_source_files", - "mcpl", "cell", "cellfrom", "cellto"), + "mcpl", "cells", "directions", "cell", "cellfrom", "cellto"), ) - if key == "surface_ids": - cv.check_type( - "surface ids for source banking", value, Iterable, Integral - ) - for surf_id in value: - cv.check_greater_than( - "surface id for source banking", surf_id, 0) - + if key in ("surface_ids", "cells"): + name = { + "surface_ids": "surface id(s) for source banking", + "cells": "Cell ID(s) for source banking", + }[key] + cv.check_type(name, value, Iterable, Integral) + for x in value: + cv.check_greater_than(name, x, 0) + elif key == "directions": + cv.check_type("directions corresponding to cells (from, to or both)", value, Iterable, str) + for direction in value: + if (direction not in ["from", "to", "both"]): + msg = "allowed values for direction: 'from', 'to', 'both' " + raise ValueError(msg) elif key == "mcpl": cv.check_type("write to an MCPL-format file", value, bool) elif key in ("max_particles", "max_source_files", "cell", "cellfrom", "cellto"): @@ -1538,19 +1552,29 @@ def _create_surf_source_read_subelement(self, root): def _create_surf_source_write_subelement(self, root): if self._surf_source_write: element = ET.SubElement(root, "surf_source_write") + if "surface_ids" in self._surf_source_write: subelement = ET.SubElement(element, "surface_ids") subelement.text = " ".join( str(x) for x in self._surf_source_write["surface_ids"] ) + if "mcpl" in self._surf_source_write: subelement = ET.SubElement(element, "mcpl") subelement.text = str(self._surf_source_write["mcpl"]).lower() + for key in ("max_particles", "max_source_files", "cell", "cellfrom", "cellto"): if key in self._surf_source_write: subelement = ET.SubElement(element, key) subelement.text = str(self._surf_source_write[key]) + for key in ("cells", "directions"): + if key in self._surf_source_write: + subelement = ET.SubElement(element, key) + subelement.text = " ".join( + str(x) for x in self._surf_source_write[key] + ) + def _create_collision_track_subelement(self, root): if self._collision_track: element = ET.SubElement(root, "collision_track") @@ -2047,9 +2071,11 @@ def _surf_source_write_from_xml_element(self, root): elem = root.find('surf_source_write') if elem is None: return - for key in ('surface_ids', 'max_particles', 'max_source_files', 'mcpl', 'cell', 'cellto', 'cellfrom'): - if key == 'surface_ids': + for key in ('surface_ids', 'max_particles', 'max_source_files', 'mcpl', 'cells', 'directions', 'cell', 'cellto', 'cellfrom'): + if key in ['surface_ids', 'cells']: value = get_elem_list(elem, key, int) + elif key == 'directions': + value = get_elem_list(elem, key, str) else: value = get_text(elem, key) if value is not None: @@ -2059,6 +2085,7 @@ def _surf_source_write_from_xml_element(self, root): value = int(value) self.surf_source_write[key] = value + def _collision_track_from_xml_element(self, root): elem = root.find('collision_track') if elem is not None: diff --git a/src/finalize.cpp b/src/finalize.cpp index 4ac6d09f321..934e815cdd7 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -129,7 +129,6 @@ int openmc_finalize() settings::source_separate = false; settings::source_write = true; settings::ssw_cell_id = C_NONE; - settings::ssw_cell_type = SSWCellType::None; settings::ssw_max_particles = 0; settings::ssw_max_files = 1; settings::survival_biasing = false; diff --git a/src/particle.cpp b/src/particle.cpp index 8d0e5b3f06e..6b38053cf11 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -889,88 +889,95 @@ void Particle::update_neutron_xs( //============================================================================== void add_surf_source_to_bank(Particle& p, const Surface& surf) { - if (simulation::current_batch <= settings::n_inactive || - simulation::surf_source_bank.full()) { + if (simulation::current_batch <= settings::n_inactive) { return; } - // If a cell/cellfrom/cellto parameter is defined - if (settings::ssw_cell_id != C_NONE) { - - // Retrieve cell index and storage type - int cell_idx = model::cell_map[settings::ssw_cell_id]; - - if (surf.bc_) { - // Leave if cellto with vacuum boundary condition - if (surf.bc_->type() == "vacuum" && - settings::ssw_cell_type == SSWCellType::To) { - return; - } + bool add_site = true; // add the site if 'cells' is not defined + // If 'cells' is defined + if (!settings::ssw_cells.empty()) { + if (surf.bc_ && surf.bc_->type() != "vacuum") { // Leave if other boundary condition than vacuum - if (surf.bc_->type() != "vacuum") { - return; - } + return; } - - // Check if the cell of interest has been exited - bool exited = false; - for (int i = 0; i < p.n_coord_last(); ++i) { - if (p.cell_last(i) == cell_idx) { - exited = true; + add_site = false; // we assume all cell-direction pairs are invalid till one + // of them passes all the tests + for (auto& cell : settings::ssw_cells) { + // Retrieve cell index and storage type + int cell_idx = model::cell_map[cell.first]; + SSWCellType direction = cell.second; + if (surf.bc_ && surf.bc_->type() == "vacuum" && + direction == SSWCellType::To) { + // skip if cellto with vacuum boundary condition + continue; } - } - // Check if the cell of interest has been entered - bool entered = false; - for (int i = 0; i < p.n_coord(); ++i) { - if (p.coord(i).cell() == cell_idx) { - entered = true; + // Check if the cell of interest has been exited + bool exited = false; + for (int i = 0; i < p.n_coord_last(); ++i) { + if (p.cell_last(i) == cell_idx) { + exited = true; + } } - } - // Vacuum boundary conditions: return if cell is not exited - if (surf.bc_) { - if (surf.bc_->type() == "vacuum" && !exited) { - return; + // Check if the cell of interest has been entered + bool entered = false; + for (int i = 0; i < p.n_coord(); ++i) { + if (p.coord(i).cell() == cell_idx) { + entered = true; + } } - } else { - // If we both enter and exit the cell of interest - if (entered && exited) { - return; - } + // Vacuum boundary conditions: return if cell is not exited + if (surf.bc_) { + if (surf.bc_->type() == "vacuum" && !exited) { + continue; + } + } else { - // If we did not enter nor exit the cell of interest - if (!entered && !exited) { - return; - } + // If we both enter and exit the cell of interest + if (entered && exited) { + continue; + } - // If cellfrom and the cell before crossing is not the cell of - // interest - if (settings::ssw_cell_type == SSWCellType::From && !exited) { - return; - } + // If we did not enter nor exit the cell of interest + if (!entered && !exited) { + continue; + } - // If cellto and the cell after crossing is not the cell of interest - if (settings::ssw_cell_type == SSWCellType::To && !entered) { - return; + // If cellfrom and the cell before crossing is not the cell of + // interest + if (direction == SSWCellType::From && !exited) { + continue; + } + + // If cellto and the cell after crossing is not the cell of interest + if (direction == SSWCellType::To && !entered) { + continue; + } } + // if a cell-direction pair survived all the checks we add the site and + // terminate the loop + add_site = true; + break; } } - SourceSite site; - site.r = p.r(); - site.u = p.u(); - site.E = p.E(); - site.time = p.time(); - site.wgt = p.wgt(); - site.delayed_group = p.delayed_group(); - site.surf_id = surf.id_; - site.particle = p.type(); - site.parent_id = p.id(); - site.progeny_id = p.n_progeny(); - int64_t idx = simulation::surf_source_bank.thread_safe_append(site); + if (add_site) { + SourceSite site; + site.r = p.r(); + site.u = p.u(); + site.E = p.E(); + site.time = p.time(); + site.wgt = p.wgt(); + site.delayed_group = p.delayed_group(); + site.surf_id = surf.id_; + site.particle = p.type(); + site.parent_id = p.id(); + site.progeny_id = p.n_progeny(); + int64_t idx = simulation::surf_source_bank.thread_safe_append(site); + } } } // namespace openmc diff --git a/src/settings.cpp b/src/settings.cpp index b6a9ab95377..023ba1b162f 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -135,7 +135,7 @@ CollisionTrackConfig collision_track_config {}; int64_t ssw_max_particles; int64_t ssw_max_files; int64_t ssw_cell_id {C_NONE}; -SSWCellType ssw_cell_type {SSWCellType::None}; +std::unordered_map ssw_cells; TemperatureMethod temperature_method {TemperatureMethod::NEAREST}; double temperature_tolerance {10.0}; double temperature_default {293.6}; @@ -925,25 +925,65 @@ void read_settings_xml(pugi::xml_node root) surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl"); } // Get cell information - if (check_for_node(node_ssw, "cell")) { - ssw_cell_id = std::stoll(get_node_value(node_ssw, "cell")); - ssw_cell_type = SSWCellType::Both; - } - if (check_for_node(node_ssw, "cellfrom")) { - if (ssw_cell_id != C_NONE) { - fatal_error( - "'cell', 'cellfrom' and 'cellto' cannot be used at the same time."); + if (check_for_node(node_ssw, "cells")) { + // raise an error if the new syntax is mixed with the old syntax + if (check_for_node(node_ssw, "cell") || + check_for_node(node_ssw, "cellfrom") || + check_for_node(node_ssw, "cellto")) { + fatal_error("'cells' cannot be used at the same time with 'cell', " + "'cellfrom' or 'cellto'."); } - ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellfrom")); - ssw_cell_type = SSWCellType::From; - } - if (check_for_node(node_ssw, "cellto")) { - if (ssw_cell_id != C_NONE) { - fatal_error( - "'cell', 'cellfrom' and 'cellto' cannot be used at the same time."); + auto ids = get_node_array(node_ssw, "cells"); + if (check_for_node(node_ssw, "directions")) { + auto directions = get_node_array(node_ssw, "directions"); + if (directions.size() != ids.size()) { + fatal_error("'directions' must have the same length as 'cells'"); + } + for (std::size_t i {0}; i < ids.size(); ++i) { + SSWCellType direction = ssw_cell_type_from_string(directions[i]); + auto [it, inserted] = ssw_cells.emplace(ids[i], direction); + // check for duplicate keys with different values + if (!inserted && it->second != direction) { + // the union of different values will always be 'Both' + it->second = SSWCellType::Both; + } + } + } else { // default behavior if 'directions' is not defined + for (std::size_t i {0}; i < ids.size(); ++i) { + ssw_cells.emplace(ids[i], SSWCellType::Both); + } + } + } else { + if (check_for_node(node_ssw, "directions")) { + fatal_error("'directions' cannot be used if 'cells' is not defined."); + } + // old syntax will be deprecated in the future + if (check_for_node(node_ssw, "cell")) { + warning("'cell' is deprecated and will be removed in the future. Use " + "'cells' and 'directions' instead."); + ssw_cell_id = std::stoll(get_node_value(node_ssw, "cell")); + ssw_cells.emplace(ssw_cell_id, SSWCellType::Both); + } + if (check_for_node(node_ssw, "cellfrom")) { + warning("'cellfrom' is deprecated and will be removed in the future. " + "Use 'cells' and 'directions' instead."); + if (ssw_cell_id != C_NONE) { + fatal_error( + "'cell', 'cellfrom' and 'cellto' cannot be used at the same time."); + } + ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellfrom")); + ssw_cells.emplace(ssw_cell_id, SSWCellType::From); + } + if (check_for_node(node_ssw, "cellto")) { + warning("'cellto' is deprecated and will be removed in the future. Use " + "'cells' and 'directions' instead."); + if (ssw_cell_id != C_NONE) { + fatal_error( + "'cell', 'cellfrom' and 'cellto' cannot be used at the same time."); + } + ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellto")); + ssw_cells.emplace(ssw_cell_id, SSWCellType::To); } - ssw_cell_id = std::stoll(get_node_value(node_ssw, "cellto")); - ssw_cell_type = SSWCellType::To; } } @@ -1277,9 +1317,21 @@ void free_memory_settings() settings::statepoint_batch.clear(); settings::sourcepoint_batch.clear(); settings::source_write_surf_id.clear(); + settings::ssw_cells.clear(); settings::res_scat_nuclides.clear(); } +SSWCellType ssw_cell_type_from_string(std::string_view s) +{ + if (s == "from") + return SSWCellType::From; + if (s == "to") + return SSWCellType::To; + if (s == "both") + return SSWCellType::Both; + throw std::invalid_argument("direction must be 'from', 'to', or 'both'"); +} + //============================================================================== // C API functions //============================================================================== diff --git a/src/surface.cpp b/src/surface.cpp index 81b756deae7..25ae5091efe 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -63,8 +63,9 @@ Surface::Surface(pugi::xml_node surf_node) { if (check_for_node(surf_node, "id")) { id_ = std::stoi(get_node_value(surf_node, "id")); - if (contains(settings::source_write_surf_id, id_) || - settings::source_write_surf_id.empty()) { + if (settings::surf_source_write && + (contains(settings::source_write_surf_id, id_) || + settings::source_write_surf_id.empty())) { surf_source_ = true; } } else { diff --git a/tests/regression_tests/surface_source_write/_visualize.py b/tests/regression_tests/surface_source_write/_visualize.py index 73340cae06f..21cb88e5f4a 100644 --- a/tests/regression_tests/surface_source_write/_visualize.py +++ b/tests/regression_tests/surface_source_write/_visualize.py @@ -10,11 +10,11 @@ # Select an option # "show": 3D visualization using matplotlib # "savefig": 2D representation using matplotlib and storing the fig under plot_2d.png - option = "show" - # option = "savefig" + #option = "show" + option = "savefig" # Select the case from its folder name - folder = "case-20" + folder = "case-24" # Reading the surface source file with h5py.File(f"{folder}/surface_source_true.h5", "r") as fp: diff --git a/tests/regression_tests/surface_source_write/case-22/inputs_true.dat b/tests/regression_tests/surface_source_write/case-22/inputs_true.dat new file mode 100644 index 00000000000..a39d905ce90 --- /dev/null +++ b/tests/regression_tests/surface_source_write/case-22/inputs_true.dat @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 60 + 5 + + + 1.5 0.5 0.5 + + + + + + 2.5 0.5 0.5 + + + + + + 0.5 0.5 1.5 + + + + + + 3.5 0.5 1.5 + + + + + + 1.5 0.5 2.5 + + + + + + 2.5 0.5 2.5 + + + + + 2 4 9 10 + 300 + 6 + to + + 1 + + diff --git a/tests/regression_tests/surface_source_write/case-22/results_true.dat b/tests/regression_tests/surface_source_write/case-22/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/surface_source_write/case-22/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-22/surface_source_true.h5 new file mode 100644 index 00000000000..27bca94959a Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-22/surface_source_true.h5 differ diff --git a/tests/regression_tests/surface_source_write/case-23/inputs_true.dat b/tests/regression_tests/surface_source_write/case-23/inputs_true.dat new file mode 100644 index 00000000000..1b36485b501 --- /dev/null +++ b/tests/regression_tests/surface_source_write/case-23/inputs_true.dat @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 60 + 5 + + + 1.5 0.5 0.5 + + + + + + 2.5 0.5 0.5 + + + + + + 0.5 0.5 1.5 + + + + + + 3.5 0.5 1.5 + + + + + + 1.5 0.5 2.5 + + + + + + 2.5 0.5 2.5 + + + + + 2 4 9 10 + 300 + 7 + to + + 1 + + diff --git a/tests/regression_tests/surface_source_write/case-23/results_true.dat b/tests/regression_tests/surface_source_write/case-23/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/surface_source_write/case-23/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-23/surface_source_true.h5 new file mode 100644 index 00000000000..432adb6a30c Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-23/surface_source_true.h5 differ diff --git a/tests/regression_tests/surface_source_write/case-24/inputs_true.dat b/tests/regression_tests/surface_source_write/case-24/inputs_true.dat new file mode 100644 index 00000000000..75fd6ab2a4e --- /dev/null +++ b/tests/regression_tests/surface_source_write/case-24/inputs_true.dat @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 60 + 5 + + + 1.5 0.5 0.5 + + + + + + 2.5 0.5 0.5 + + + + + + 0.5 0.5 1.5 + + + + + + 3.5 0.5 1.5 + + + + + + 1.5 0.5 2.5 + + + + + + 2.5 0.5 2.5 + + + + + 2 4 9 10 + 300 + 6 7 + to to + + 1 + + diff --git a/tests/regression_tests/surface_source_write/case-24/results_true.dat b/tests/regression_tests/surface_source_write/case-24/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/surface_source_write/case-24/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-24/surface_source_true.h5 new file mode 100644 index 00000000000..5fd380aba2d Binary files /dev/null and b/tests/regression_tests/surface_source_write/case-24/surface_source_true.h5 differ diff --git a/tests/regression_tests/surface_source_write/test.py b/tests/regression_tests/surface_source_write/test.py index 094df1f8b84..55d118eec71 100644 --- a/tests/regression_tests/surface_source_write/test.py +++ b/tests/regression_tests/surface_source_write/test.py @@ -25,6 +25,7 @@ - model_2: cylindrical core in 1 box (vacuum BC), - model_3: cylindrical core in 1 box (reflective BC), - model_4: cylindrical core in 1 box (periodic BC). +- model_5: 4*1*3 array of boxes (transmission BC) Two models including DAGMC geometries are also used, based on the mesh file 'dagmc.h5m' available from tests/regression_tests/dagmc/legacy: @@ -79,6 +80,12 @@ case-20 model_4 1 No P+R Particles crossing the declared periodic surface case-21 model_4 1 cell (root universe) P+R None +case-22 model_5 Multiple cellto T particles crossing the declared + surfaces that enters a cell +case-23 model_5 Multiple cellto T particles crossing the declared + surfaces that enters a cell +case-24 model_5 Multiple cellto (multiple) T particles crossing the declared + surfaces that enters multiple cells ======== ======= ========= ========================= ===== =================================== *: BC stands for Boundary Conditions, T for Transmission, R for Reflective, and V for Vacuum. @@ -603,6 +610,80 @@ def model_4(): return model +@pytest.fixture +def model_5(): + """4*1*3 array of boxes""" + openmc.reset_auto_ids() + model = openmc.Model() + + # ============================================================================= + # Materials + # ============================================================================= + + # ============================================================================= + # Geometry + # ============================================================================= + + nx = 4 + ny = 1 + nz = 3 + + x_planes = [None] * (nx+1) + y_planes = [None] * (ny+1) + z_planes = [None] * (nz+1) + + for i in range(nx+1): + x_planes[i] = openmc.XPlane(x0=i) + for i in range(ny+1): + y_planes[i] = openmc.YPlane(y0=i) + for i in range(nz+1): + z_planes[i] = openmc.ZPlane(z0=i) + + for planes in (x_planes, y_planes, z_planes): + for i in (0,-1): + planes[i].boundary_type = 'vacuum' + + cells = [[None]*nz for _ in range(nx)] + for j in range(nz): + for i in range(nx): + cells[i][j] = openmc.Cell(region = +x_planes[i] & -x_planes[i+1] & +y_planes[0] & -y_planes[-1] & +z_planes[j] & -z_planes[j+1]) + + cells_1D = [cells[i][j] for j in range(len(cells[0])) for i in range(len(cells))] + root = openmc.Universe(cells=cells_1D) + model.geometry = openmc.Geometry(root) + + # ============================================================================= + # Settings + # ============================================================================= + + model.settings = openmc.Settings() + model.settings.run_mode = 'fixed source' + model.settings.particles = 60 + model.settings.batches = 5 + model.settings.seed = 1 + + point_21 = openmc.stats.Point((1.5,0.5,0.5)) + point_31 = openmc.stats.Point((2.5,0.5,0.5)) + point_12 = openmc.stats.Point((0.5,0.5,1.5)) + point_42 = openmc.stats.Point((3.5,0.5,1.5)) + point_23 = openmc.stats.Point((1.5,0.5,2.5)) + point_33 = openmc.stats.Point((2.5,0.5,2.5)) + + x_pos = openmc.stats.Monodirectional((1,0,0)) + x_neg = openmc.stats.Monodirectional((-1,0,0)) + z_pos = openmc.stats.Monodirectional((0,0,1)) + z_neg = openmc.stats.Monodirectional((0,0,-1)) + + source_1 = openmc.IndependentSource(space=point_21, angle=z_pos, strength=1.0) + source_2 = openmc.IndependentSource(space=point_31, angle=z_pos, strength=1.0) + source_3 = openmc.IndependentSource(space=point_12, angle=x_pos, strength=1.0) + source_4 = openmc.IndependentSource(space=point_42, angle=x_neg, strength=1.0) + source_5 = openmc.IndependentSource(space=point_23, angle=z_neg, strength=1.0) + source_6 = openmc.IndependentSource(space=point_33, angle=z_neg, strength=1.0) + model.settings.source = [source_1, source_2, source_3, source_4, source_5, source_6] + + return model + def return_surface_source_data(filepath): """Read a surface source file and return a sorted array composed @@ -817,6 +898,21 @@ def _cleanup(self): "model_4", {"max_particles": 300, "surface_ids": [4], "cell": 3}, ), + ( + "case-22", + "model_5", + {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [6], "directions": ["to"]}, + ), + ( + "case-23", + "model_5", + {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [7], "directions": ["to"]}, + ), + ( + "case-24", + "model_5", + {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [6, 7], "directions": ["to", "to"]}, + ), ], ) def test_surface_source_cell_history_based( @@ -1123,3 +1219,83 @@ def test_surface_source_cell_dagmc( "statepoint.5.h5", model=model, workdir=folder ) harness.main() + +def test_surface_source_multiple_cells(model_5, single_thread, single_process): + """Test that the number of particles entering two cells equal the sum of the number of particles entering each individual cell""" + assert os.environ["OMP_NUM_THREADS"] == "1" + assert config["mpi_np"] == "1" + + def run_and_count(model, folder, parameter): + model.settings.surf_source_write = parameter + harness = SurfaceSourceWriteTestHarness("statepoint.5.h5", model=model, workdir=folder) + + base_dir = os.getcwd() + try: + os.chdir(folder) + harness._build_inputs() + inputs = harness._get_inputs() + harness._write_inputs(inputs) + harness._compare_inputs() + harness._run_openmc() + harness._test_output_created() + return len(return_surface_source_data("surface_source.h5")) + finally: + harness._cleanup() + os.chdir(base_dir) + + p1 = {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [6], "directions": ["to"]} + p2 = {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [7], "directions": ["to"]} + p_sum = {"max_particles": 300, "surface_ids": [2, 4, 9, 10], "cells": [6, 7], "directions": ["to", "to"]} + + n1 = run_and_count(model_5, "case-22", p1) + n2 = run_and_count(model_5, "case-23", p2) + n_sum = run_and_count(model_5, "case-24", p_sum) + + assert n_sum == n1 + n2 + +def test_duplicate_cells(tmp_path, model_5, single_thread, single_process): + """Test the equivalence of inputs having duplicate cells""" + assert os.environ["OMP_NUM_THREADS"] == "1" + assert config["mpi_np"] == "1" + + params = [ + {"max_particles": 700, "surface_ids": [2, 4, 9, 10], + "cells": [6, 7, 7], "directions": ["both", "from", "to"]}, + + {"max_particles": 700, "surface_ids": [2, 4, 9, 10], + "cells": [6, 7, 7], "directions": ["both", "both", "to"]}, + + {"max_particles": 700, "surface_ids": [2, 4, 9, 10], + "cells": [6, 7, 7], "directions": ["both", "to", "both"]}, + + {"max_particles": 700, "surface_ids": [2, 4, 9, 10], + "cells": [6, 7, 7, 7], "directions": ["both", "both", "to", "from"]}, + ] + + def run_and_read(parameter, subdir): + run_dir = tmp_path / subdir + run_dir.mkdir() + + model = model_5 + model.settings.surf_source_write = parameter + + harness = SurfaceSourceWriteTestHarness("statepoint.5.h5", model=model, workdir=str(run_dir)) + + base = os.getcwd() + try: + os.chdir(run_dir) + harness._build_inputs() + inputs = harness._get_inputs() + harness._write_inputs(inputs) + harness._run_openmc() + harness._test_output_created() + return return_surface_source_data("surface_source.h5") + finally: + harness._cleanup() + os.chdir(base) + + ref = run_and_read(params[0], "run0") + for k, p in enumerate(params[1:], start=1): + out = run_and_read(p, f"run{k}") + assert out.shape == ref.shape + assert np.array_equal(out, ref) \ No newline at end of file