From 7fc5c7e5479429815b8f3f7ff350715c730a5e46 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 27 Sep 2025 20:51:09 -0500 Subject: [PATCH 01/11] use new API for adaptive mesh --- include/openmc/mesh.h | 19 +++++++++++++++++-- src/mesh.cpp | 42 +++++++++++++++++++++++++++++++++++++++++- 2 files changed, 58 insertions(+), 3 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 2936d3ab1bd..a5f942548a4 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -429,7 +429,6 @@ class PeriodicStructuredMesh : public StructuredMesh { public: PeriodicStructuredMesh() = default; PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; - PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {}; Position local_coords(const Position& r) const override { @@ -788,7 +787,6 @@ class MOABMesh : public UnstructuredMesh { // Constructors MOABMesh() = default; MOABMesh(pugi::xml_node); - MOABMesh(hid_t group); MOABMesh(const std::string& filename, double length_multiplier = 1.0); MOABMesh(std::shared_ptr external_mbi); @@ -1062,6 +1060,9 @@ class AdaptiveLibMesh : public LibMesh { int get_bin(Position r) const override; + //! setter for mesh tally amalgamtion + void set_mesh_tally_amalgamation(std::string cluster_element_integer_name); + protected: // Overridden methods int get_bin_from_element(const libMesh::Elem* elem) const override; @@ -1080,6 +1081,20 @@ class AdaptiveLibMesh : public LibMesh { //!< elements std::vector elem_to_bin_map_; //!< mapping dof indices to bin indices for //!< active elements + + bool amalgamation_ = false; //!< whether we are doing mesh and tally + //!< amalgamation by default it's turned off. + + int clustering_element_integer_index_ = -1; //!< extra element integer index for + // element clustering + + /*create a hash map where every element in a cluster would map to the first + * element of in that cluster if the element isn't part of a cluster then it + * will point to it self + */ + std::unordered_map + clustering_element_mapping_; }; #endif diff --git a/src/mesh.cpp b/src/mesh.cpp index 6ee15b9c18b..e506bdbc610 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3763,6 +3763,44 @@ AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, } } +void AdaptiveLibMesh::set_mesh_tally_amalgamation( + std::string cluster_element_integer_name) +{ + + clustering_element_integer_index_ = + m_->has_elem_integer(cluster_element_integer_name) + ? m_->get_elem_integer_index(cluster_element_integer_name) + : -1; + amalgamation_ = (clustering_element_integer_index_ != -1); + + if (amalgamation_) { + // reseve the hash map for cluster elements + clustering_element_mapping_.reserve(m_->n_active_elem()); + + // adding clustering map + for (auto it = m_->active_elements_begin(); it != m_->active_elements_end(); + it++) { + + auto elem = *it; + auto cluster_elem = elem; + unsigned int cluster_id = + elem->get_extra_integer(clustering_element_integer_index_); + + if (cluster_id != -1) { + auto first_element_in_a_cluster = m_->elem_ptr(cluster_id); + + if (first_element_in_a_cluster and first_element_in_a_cluster->active()) + cluster_elem = first_element_in_a_cluster; + } + clustering_element_mapping_.insert(std::make_pair(elem, cluster_elem)); + } + } + else{ + fatal_error(fmt::format("No extra element integer named: {} found in the " + "mesh", cluster_element_integer_name)); + } +} + int AdaptiveLibMesh::n_bins() const { return num_active_; @@ -3813,7 +3851,9 @@ int AdaptiveLibMesh::get_bin(Position r) const int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const { - int bin = elem_to_bin_map_[elem->id()]; + auto tally_elem = amalgamation_ ? clustering_element_mapping_.at(elem) : elem; + int bin = adaptive_ ? elem_to_bin_map_[tally_elem->id()] + : tally_elem->id() - first_element_id_; if (bin >= n_bins() || bin < 0) { fatal_error(fmt::format("Invalid bin: {}", bin)); } From 4ff4d4cdfd7aecdb4c6667c78855b906915c9c38 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 27 Sep 2025 21:43:06 -0500 Subject: [PATCH 02/11] I don't need adaptive_ check for an Adaptive class --- include/openmc/mesh.h | 5 +++-- src/mesh.cpp | 19 ++++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index a5f942548a4..a605d59812b 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -1085,8 +1085,9 @@ class AdaptiveLibMesh : public LibMesh { bool amalgamation_ = false; //!< whether we are doing mesh and tally //!< amalgamation by default it's turned off. - int clustering_element_integer_index_ = -1; //!< extra element integer index for - // element clustering + int clustering_element_integer_index_ = + -1; //!< extra element integer index for + // element clustering /*create a hash map where every element in a cluster would map to the first * element of in that cluster if the element isn't part of a cluster then it diff --git a/src/mesh.cpp b/src/mesh.cpp index e506bdbc610..a3783dbdc89 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3766,20 +3766,20 @@ AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, void AdaptiveLibMesh::set_mesh_tally_amalgamation( std::string cluster_element_integer_name) { - + clustering_element_integer_index_ = m_->has_elem_integer(cluster_element_integer_name) ? m_->get_elem_integer_index(cluster_element_integer_name) : -1; amalgamation_ = (clustering_element_integer_index_ != -1); - + if (amalgamation_) { // reseve the hash map for cluster elements clustering_element_mapping_.reserve(m_->n_active_elem()); // adding clustering map for (auto it = m_->active_elements_begin(); it != m_->active_elements_end(); - it++) { + it++) { auto elem = *it; auto cluster_elem = elem; @@ -3794,10 +3794,10 @@ void AdaptiveLibMesh::set_mesh_tally_amalgamation( } clustering_element_mapping_.insert(std::make_pair(elem, cluster_elem)); } - } - else{ + } else { fatal_error(fmt::format("No extra element integer named: {} found in the " - "mesh", cluster_element_integer_name)); + "mesh", + cluster_element_integer_name)); } } @@ -3852,9 +3852,10 @@ int AdaptiveLibMesh::get_bin(Position r) const int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const { auto tally_elem = amalgamation_ ? clustering_element_mapping_.at(elem) : elem; - int bin = adaptive_ ? elem_to_bin_map_[tally_elem->id()] - : tally_elem->id() - first_element_id_; - if (bin >= n_bins() || bin < 0) { + int bin = elem_to_bin_map_[tally_elem->id()]; + + if (bin >= n_bins() || bin < 0) + { fatal_error(fmt::format("Invalid bin: {}", bin)); } return bin; From 5c167568f059f2611dd3ed2371690d14463a1349 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sun, 28 Sep 2025 08:59:15 -0500 Subject: [PATCH 03/11] apply clang format --- src/mesh.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index a3783dbdc89..c909993b76e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3854,8 +3854,7 @@ int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const auto tally_elem = amalgamation_ ? clustering_element_mapping_.at(elem) : elem; int bin = elem_to_bin_map_[tally_elem->id()]; - if (bin >= n_bins() || bin < 0) - { + if (bin >= n_bins() || bin < 0) { fatal_error(fmt::format("Invalid bin: {}", bin)); } return bin; From 752e80b2e8eaa1abe313c2309343af63ee3e4c00 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 17 Jan 2026 11:10:04 -0600 Subject: [PATCH 04/11] applying clang format --- include/openmc/mesh.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index a605d59812b..2ab8c412871 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -429,6 +429,7 @@ class PeriodicStructuredMesh : public StructuredMesh { public: PeriodicStructuredMesh() = default; PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {}; + PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {}; Position local_coords(const Position& r) const override { @@ -787,6 +788,7 @@ class MOABMesh : public UnstructuredMesh { // Constructors MOABMesh() = default; MOABMesh(pugi::xml_node); + MOABMesh(hid_t group); MOABMesh(const std::string& filename, double length_multiplier = 1.0); MOABMesh(std::shared_ptr external_mbi); From 5d2c6973d3fce6415f88e585fa15fbc69f8758ec Mon Sep 17 00:00:00 2001 From: "Travis L." Date: Mon, 27 Apr 2026 08:40:44 -0600 Subject: [PATCH 05/11] Clean up MCPL references (#3927) --- docs/source/usersguide/install.rst | 10 +++++----- src/output.cpp | 5 ----- 2 files changed, 5 insertions(+), 10 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 1ef2f574be8..2a0d301d4bc 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -331,11 +331,11 @@ Prerequisites This option allows OpenMC to read and write MCPL (Monte Carlo Particle Lists) files instead of .h5 files for sources (external source - distribution, k-eigenvalue source distribution, and surface sources). To - turn this option on in the CMake configuration step, add the following - option:: - - cmake -DOPENMC_USE_MCPL=on .. + distribution, k-eigenvalue source distribution, and surface sources). + OpenMC does not need any particular build option to use this, but MCPL + must be installed on the system in order to do so. Refer to the + `MCPL documentation `_ + for instructions on how to accomplish this. * NCrystal_ library for defining materials with enhanced thermal neutron transport diff --git a/src/output.cpp b/src/output.cpp index c66c313dce5..f8c0f2a97a7 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -321,7 +321,6 @@ void print_build_info() std::string png(n); std::string profiling(n); std::string coverage(n); - std::string mcpl(n); std::string uwuw(n); std::string strict_fp(n); @@ -337,9 +336,6 @@ void print_build_info() #ifdef OPENMC_LIBMESH_ENABLED libmesh = y; #endif -#ifdef OPENMC_MCPL - mcpl = y; -#endif #ifdef USE_LIBPNG png = y; #endif @@ -369,7 +365,6 @@ void print_build_info() fmt::print("PNG support: {}\n", png); fmt::print("DAGMC support: {}\n", dagmc); fmt::print("libMesh support: {}\n", libmesh); - fmt::print("MCPL support: {}\n", mcpl); fmt::print("Coverage testing: {}\n", coverage); fmt::print("Profiling flags: {}\n", profiling); fmt::print("UWUW support: {}\n", uwuw); From 8c6f4e21caebd12991a307ad6e8a688331ca0273 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Apr 2026 09:02:08 -0500 Subject: [PATCH 06/11] Support multiple meshes in R2S calculations (#3860) --- docs/source/usersguide/decay_sources.rst | 19 ++ openmc/deplete/microxs.py | 145 ++++++++------ openmc/deplete/r2s.py | 242 +++++++++++++---------- tests/unit_tests/test_deplete_microxs.py | 6 +- tests/unit_tests/test_r2s.py | 75 ++++++- 5 files changed, 316 insertions(+), 171 deletions(-) diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index 398680e7464..21981fdaa82 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -190,6 +190,25 @@ we would run:: r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000}) +It is also possible to use multiple meshes by passing a list of meshes instead +of a single mesh. This can be useful, for example, when different regions of the +model require different mesh resolutions. The meshes are assumed to be +**non-overlapping**; each element--material combination across all meshes is +treated as an independent activation region, and all meshes are handled in a +single neutron transport solve. For example:: + + # Fine mesh near the activation target + mesh_fine = openmc.RegularMesh() + mesh_fine.dimension = (10, 10, 10) + ... + + # Coarse mesh for the surrounding region + mesh_coarse = openmc.RegularMesh() + mesh_coarse.dimension = (5, 5, 5) + ... + + r2s = openmc.deplete.R2SManager(model, [mesh_fine, mesh_coarse]) + Direct 1-Step (D1S) Calculations ================================ diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index e6e2dbce24c..687cf646f29 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -36,7 +36,8 @@ Sequence[openmc.Cell], Sequence[openmc.Universe], openmc.MeshBase, - openmc.Filter + openmc.Filter, + Sequence[openmc.Filter] ] @@ -69,8 +70,12 @@ def get_microxs_and_flux( ---------- model : openmc.Model OpenMC model object. Must contain geometry, materials, and settings. - domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter + domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter, or list of openmc.Filter Domains in which to tally reaction rates, or a spatial tally filter. + A list of filters can be provided to create one set of tallies per + filter (e.g., one :class:`~openmc.MeshMaterialFilter` per mesh) that + are all evaluated in a single transport solve. Results are + concatenated across all filters in order. nuclides : list of str Nuclides to get cross sections for. If not specified, all burnable nuclides from the depletion chain file are used. @@ -142,26 +147,24 @@ def get_microxs_and_flux( else: energy_filter = openmc.EnergyFilter(energies) + # Build list of domain filters if isinstance(domains, openmc.Filter): - domain_filter = domains + domain_filters = [domains] elif isinstance(domains, openmc.MeshBase): - domain_filter = openmc.MeshFilter(domains) + domain_filters = [openmc.MeshFilter(domains)] + elif isinstance(domains, Sequence) and len(domains) > 0 and \ + isinstance(domains[0], openmc.Filter): + domain_filters = list(domains) elif isinstance(domains[0], openmc.Material): - domain_filter = openmc.MaterialFilter(domains) + domain_filters = [openmc.MaterialFilter(domains)] elif isinstance(domains[0], openmc.Cell): - domain_filter = openmc.CellFilter(domains) + domain_filters = [openmc.CellFilter(domains)] elif isinstance(domains[0], openmc.Universe): - domain_filter = openmc.UniverseFilter(domains) + domain_filters = [openmc.UniverseFilter(domains)] else: raise ValueError(f"Unsupported domain type: {type(domains[0])}") - flux_tally = openmc.Tally(name='MicroXS flux') - flux_tally.filters = [domain_filter, energy_filter] - flux_tally.scores = ['flux'] - model.tallies = [flux_tally] - - # Prepare reaction-rate tally for 'direct' or subset for 'flux' with opts - rr_tally = None + # Prepare reaction-rate nuclides/reactions rr_nuclides: list[str] = [] rr_reactions: list[str] = [] if reaction_rate_mode == 'direct': @@ -177,20 +180,33 @@ def get_microxs_and_flux( if rr_reactions: rr_reactions = [r for r in rr_reactions if r in set(reactions)] - # Only construct tally if both lists are non-empty - if rr_nuclides and rr_reactions: - rr_tally = openmc.Tally(name='MicroXS RR') - # Use 1-group energy filter for RR in flux mode - if reaction_rate_mode == 'flux': - rr_energy_filter = openmc.EnergyFilter( - [energy_filter.values[0], energy_filter.values[-1]]) - else: - rr_energy_filter = energy_filter - rr_tally.filters = [domain_filter, rr_energy_filter] - rr_tally.nuclides = rr_nuclides - rr_tally.multiply_density = False - rr_tally.scores = rr_reactions - model.tallies.append(rr_tally) + # Use 1-group energy filter for RR in flux mode + has_rr = bool(rr_nuclides and rr_reactions) + if has_rr and reaction_rate_mode == 'flux': + rr_energy_filter = openmc.EnergyFilter( + [energy_filter.values[0], energy_filter.values[-1]]) + else: + rr_energy_filter = energy_filter + + # Create one flux tally (and optionally one RR tally) per domain filter. + flux_tallies = [] + rr_tallies = [] + model.tallies = [] + for i, domain_filter in enumerate(domain_filters): + flux_tally = openmc.Tally(name=f'MicroXS flux {i}') + flux_tally.filters = [domain_filter, energy_filter] + flux_tally.scores = ['flux'] + model.tallies.append(flux_tally) + flux_tallies.append(flux_tally) + + if has_rr: + rr_tally = openmc.Tally(name=f'MicroXS RR {i}') + rr_tally.filters = [domain_filter, rr_energy_filter] + rr_tally.nuclides = rr_nuclides + rr_tally.multiply_density = False + rr_tally.scores = rr_reactions + model.tallies.append(rr_tally) + rr_tallies.append(rr_tally) if openmc.lib.is_initialized: openmc.lib.finalize() @@ -227,40 +243,41 @@ def get_microxs_and_flux( # Read in tally results (on all ranks) with StatePoint(statepoint_path) as sp: - if rr_tally is not None: - rr_tally = sp.tallies[rr_tally.id] - rr_tally._read_results() - flux_tally = sp.tallies[flux_tally.id] - flux_tally._read_results() - - # Get flux values and make energy groups last dimension - flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1) - flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups) - - # Create list where each item corresponds to one domain - fluxes = list(flux.squeeze((1, 2))) - - # If we built a reaction-rate tally, compute microscopic cross sections - if rr_tally is not None: - # Get reaction rates - reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions) - - # Make energy groups last dimension - reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups) - - # If RR is 1-group, sum flux over groups - if reaction_rate_mode == "flux": - flux = flux.sum(axis=-1, keepdims=True) # (domains, 1, 1, 1) - - # Divide RR by flux to get microscopic cross sections. The indexing - # ensures that only non-zero flux values are used, and broadcasting is - # applied to align the shapes of reaction_rates and flux for division. - xs = np.zeros_like(reaction_rates) # (domains, nuclides, reactions, groups) - d, _, _, g = np.nonzero(flux) - xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g] - - # Create lists where each item corresponds to one domain - direct_micros = [MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs] + for i in range(len(flux_tallies)): + flux_tallies[i] = sp.tallies[flux_tallies[i].id] + flux_tallies[i]._read_results() + if rr_tallies: + rr_tallies[i] = sp.tallies[rr_tallies[i].id] + rr_tallies[i]._read_results() + + # Concatenate results across all domain filters + fluxes = [] + all_flux_arrays = [] + for flux_tally in flux_tallies: + # Get flux values and make energy groups last dimension + flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1) + flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups) + all_flux_arrays.append(flux) + fluxes.extend(flux.squeeze((1, 2))) + + # If we built reaction-rate tallies, compute microscopic cross sections + if rr_tallies: + direct_micros = [] + for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies): + flux = flux_arr + # Get reaction rates and make energy groups last dimension + reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions) + reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups) + + # If RR is 1-group, sum flux over groups + if reaction_rate_mode == "flux": + flux = flux.sum(axis=-1, keepdims=True) + + xs = np.zeros_like(reaction_rates) + d, _, _, g = np.nonzero(flux) + xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g] + direct_micros.extend( + MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs) # If using flux mode, compute flux-collapsed microscopic XS if reaction_rate_mode == 'flux': @@ -273,9 +290,9 @@ def get_microxs_and_flux( ) for flux_i in fluxes] # Decide which micros to use and merge if needed - if reaction_rate_mode == 'flux' and rr_tally is not None: + if reaction_rate_mode == 'flux' and rr_tallies: micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)] - elif rr_tally is not None: + elif rr_tallies: micros = direct_micros else: micros = flux_micros diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 57bbe437ffb..10d7fdd502a 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -13,11 +13,11 @@ from ..checkvalue import PathLike from ..mpi import comm from openmc.lib import TemporarySession -from openmc.utility_funcs import change_directory def get_activation_materials( - model: openmc.Model, mmv: openmc.MeshMaterialVolumes + model: openmc.Model, + mmv_list: list[openmc.MeshMaterialVolumes] ) -> openmc.Materials: """Get a list of activation materials for each mesh element/material. @@ -31,35 +31,35 @@ def get_activation_materials( ---------- model : openmc.Model The full model containing the geometry and materials. - mmv : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element. + mmv_list : list of openmc.MeshMaterialVolumes + List of mesh material volumes objects, one per mesh, containing the + materials and their volumes for each mesh element. Returns ------- openmc.Materials A list of materials, each corresponding to a unique mesh element and - material combination. + material combination across all meshes. """ - # Get the material ID, volume, and element index for each element-material - # combination - mat_ids = mmv._materials[mmv._materials > -1] - volumes = mmv._volumes[mmv._materials > -1] - elems, _ = np.where(mmv._materials > -1) - # Get all materials in the model material_dict = model._get_all_materials() # Create a new activation material for each element-material combination + # across all meshes materials = openmc.Materials() - for elem, mat_id, vol in zip(elems, mat_ids, volumes): - mat = material_dict[mat_id] - new_mat = mat.clone() - new_mat.depletable = True - new_mat.name = f'Element {elem}, Material {mat_id}' - new_mat.volume = vol - materials.append(new_mat) + for mesh_idx, mmv in enumerate(mmv_list): + mat_ids = mmv._materials[mmv._materials > -1] + volumes = mmv._volumes[mmv._materials > -1] + elems, _ = np.where(mmv._materials > -1) + + for elem, mat_id, vol in zip(elems, mat_ids, volumes): + mat = material_dict[mat_id] + new_mat = mat.clone() + new_mat.depletable = True + new_mat.name = f'Mesh {mesh_idx}, Element {elem}, Material {mat_id}' + new_mat.volume = vol + materials.append(new_mat) return materials @@ -70,7 +70,9 @@ class R2SManager: This class is responsible for managing the materials and sources needed for mesh-based or cell-based R2S calculations. It provides methods to get activation materials and decay photon sources based on the mesh/cells and - materials in the OpenMC model. + materials in the OpenMC model. Multiple meshes can be specified as domains, + in which case each element--material combination of each mesh is treated as + an activation region (meshes are assumed to be non-overlapping). This class supports the use of a different models for the neutron and photon transport calculation. However, for cell-based calculations, it assumes that @@ -83,17 +85,20 @@ class R2SManager: ---------- neutron_model : openmc.Model The OpenMC model to use for neutron transport. - domains : openmc.MeshBase or Sequence[openmc.Cell] - The mesh or a sequence of cells that represent the spatial units over - which the R2S calculation will be performed. + domains : openmc.MeshBase or Sequence[openmc.MeshBase] or Sequence[openmc.Cell] + The mesh(es) or a sequence of cells that represent the spatial units + over which the R2S calculation will be performed. When a single + :class:`~openmc.MeshBase` or a sequence of meshes is given, each + element--material combination across all meshes is treated as an + activation region. photon_model : openmc.Model, optional The OpenMC model to use for photon transport calculations. If None, a shallow copy of the neutron_model will be created and used. Attributes ---------- - domains : openmc.MeshBase or Sequence[openmc.Cell] - The mesh or a sequence of cells that represent the spatial units over + domains : list of openmc.MeshBase or Sequence[openmc.Cell] + The meshes or a sequence of cells that represent the spatial units over which the R2S calculation will be performed. neutron_model : openmc.Model The OpenMC model used for neutron transport. @@ -101,7 +106,7 @@ class R2SManager: The OpenMC model used for photon transport calculations. method : {'mesh-based', 'cell-based'} Indicates whether the R2S calculation uses mesh elements ('mesh-based') - as the spatial discetization or a list of a cells ('cell-based'). + as the spatial discretization or a list of cells ('cell-based'). results : dict A dictionary that stores results from the R2S calculation. @@ -109,7 +114,7 @@ class R2SManager: def __init__( self, neutron_model: openmc.Model, - domains: openmc.MeshBase | Sequence[openmc.Cell], + domains: openmc.MeshBase | Sequence[openmc.MeshBase] | Sequence[openmc.Cell], photon_model: openmc.Model | None = None, ): self.neutron_model = neutron_model @@ -126,9 +131,14 @@ def __init__( self.photon_model = photon_model if isinstance(domains, openmc.MeshBase): self.method = 'mesh-based' + self.domains = [domains] + elif isinstance(domains, Sequence) and len(domains) > 0 and \ + isinstance(domains[0], openmc.MeshBase): + self.method = 'mesh-based' + self.domains = list(domains) else: self.method = 'cell-based' - self.domains = domains + self.domains = list(domains) self.results = {} def run( @@ -243,11 +253,13 @@ def step1_neutron_transport( ): """Run the neutron transport step. - This step computes the material volume fractions on the mesh, creates a - mesh-material filter, and retrieves the fluxes and microscopic cross - sections for each mesh/material combination. This step will populate the - 'fluxes' and 'micros' keys in the results dictionary. For a mesh-based - calculation, it will also populate the 'mesh_material_volumes' key. + This step computes the material volume fractions on each mesh, creates + mesh-material filters, and retrieves the fluxes and microscopic cross + sections for each mesh/material combination via a single transport + solve. This step will populate the 'fluxes' and 'micros' keys in the + results dictionary. For a mesh-based calculation, it will also populate + the 'mesh_material_volumes' key (a list of + :class:`~openmc.MeshMaterialVolumes`, one per mesh). Parameters ---------- @@ -266,19 +278,28 @@ def step1_neutron_transport( output_dir.mkdir(parents=True, exist_ok=True) if self.method == 'mesh-based': - # Compute material volume fractions on the mesh + # Compute material volume fractions on each mesh if mat_vol_kwargs is None: mat_vol_kwargs = {} mat_vol_kwargs.setdefault('bounding_boxes', True) - self.results['mesh_material_volumes'] = mmv = comm.bcast( - self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs)) - # Save results to file - if comm.rank == 0: - mmv.save(output_dir / 'mesh_material_volumes.npz') + mmv_list = [] + domain_filters = [] + for i, mesh in enumerate(self.domains): + mmv = comm.bcast( + mesh.material_volumes(self.neutron_model, **mat_vol_kwargs)) + mmv_list.append(mmv) + + # Save results to file + if comm.rank == 0: + mmv.save(output_dir / f'mesh_material_volumes_{i}.npz') + + # Create mesh-material filter for this mesh + domain_filters.append( + openmc.MeshMaterialFilter.from_volumes(mesh, mmv)) - # Create mesh-material filter based on what combos were found - domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv) + self.results['mesh_material_volumes'] = mmv_list + domains = domain_filters else: domains: Sequence[openmc.Cell] = self.domains @@ -357,8 +378,9 @@ def step2_activation( if self.method == 'mesh-based': # Get unique material for each (mesh, material) combination - mmv = self.results['mesh_material_volumes'] - self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv) + mmv_list = self.results['mesh_material_volumes'] + self.results['activation_materials'] = get_activation_materials( + self.neutron_model, mmv_list) else: # Create unique material for each cell activation_mats = openmc.Materials() @@ -468,12 +490,20 @@ def step3_photon_transport( # photon model if it is different from the neutron model to account for # potential material changes if self.method == 'mesh-based' and different_photon_model: - self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast( - self.domains.material_volumes(self.photon_model, **mat_vol_kwargs)) + if mat_vol_kwargs is None: + mat_vol_kwargs = {} + photon_mmv_list = [] + for i, mesh in enumerate(self.domains): + photon_mmv = comm.bcast( + mesh.material_volumes(self.photon_model, **mat_vol_kwargs)) + photon_mmv_list.append(photon_mmv) + + # Save photon MMV results to file + if comm.rank == 0: + photon_mmv.save( + output_dir / f'mesh_material_volumes_{i}.npz') - # Save photon MMV results to file - if comm.rank == 0: - photon_mmv.save(output_dir / 'mesh_material_volumes.npz') + self.results['mesh_material_volumes_photon'] = photon_mmv_list if comm.rank == 0: tally_ids = [tally.id for tally in self.photon_model.tallies] @@ -543,7 +573,7 @@ def get_decay_photon_source_mesh( ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. - For each mesh element-material combination, an + For each mesh element-material combination across all meshes, an :class:`~openmc.IndependentSource` is created with a :class:`~openmc.stats.Box` spatial distribution based on the bounding box of the material within the mesh element. A material constraint is @@ -575,52 +605,56 @@ def get_decay_photon_source_mesh( index_mat = 0 # Get various results from previous steps - mat_vols = self.results['mesh_material_volumes'] + mmv_list = self.results['mesh_material_volumes'] materials = self.results['activation_materials'] results = self.results['depletion_results'] - photon_mat_vols = self.results.get('mesh_material_volumes_photon') - - # Total number of mesh elements - n_elements = mat_vols.num_elements - - for index_elem in range(n_elements): - # Determine which materials exist in the photon model for this element - if photon_mat_vols is not None: - photon_materials = { - mat_id - for mat_id, _ in photon_mat_vols.by_element(index_elem) - if mat_id is not None - } - - for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): - # Skip void volume - if mat_id is None: - continue - - # Skip if this material doesn't exist in photon model - if photon_mat_vols is not None and mat_id not in photon_materials: - index_mat += 1 - continue - - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source - energy = activated_mat.get_decay_photon_energy() - if energy is not None: - strength = energy.integral() - space = openmc.stats.Box(*bbox) - sources.append(openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [mat_dict[mat_id]]} - )) + photon_mmv_list = self.results.get('mesh_material_volumes_photon') + + for mesh_idx, mat_vols in enumerate(mmv_list): + photon_mat_vols = photon_mmv_list[mesh_idx] \ + if photon_mmv_list is not None else None + + # Total number of mesh elements for this mesh + n_elements = mat_vols.num_elements + + for index_elem in range(n_elements): + # Determine which materials exist in the photon model for this element + if photon_mat_vols is not None: + photon_materials = { + mat_id + for mat_id, _ in photon_mat_vols.by_element(index_elem) + if mat_id is not None + } + + for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): + # Skip void volume + if mat_id is None: + continue + + # Skip if this material doesn't exist in photon model + if photon_mat_vols is not None and mat_id not in photon_materials: + index_mat += 1 + continue + + # Get activated material composition + original_mat = materials[index_mat] + activated_mat = results[time_index].get_material(str(original_mat.id)) - # Increment index of activated material - index_mat += 1 + # Create decay photon source + energy = activated_mat.get_decay_photon_energy() + if energy is not None: + strength = energy.integral() + space = openmc.stats.Box(*bbox) + sources.append(openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=strength, + constraints={'domains': [mat_dict[mat_id]]} + )) + + # Increment index of activated material + index_mat += 1 return sources @@ -638,10 +672,13 @@ def load_results(self, path: PathLike): # Load neutron transport results neutron_dir = path / 'neutron_transport' if self.method == 'mesh-based': - mmv_file = neutron_dir / 'mesh_material_volumes.npz' - if mmv_file.exists(): - self.results['mesh_material_volumes'] = \ - openmc.MeshMaterialVolumes.from_npz(mmv_file) + mmv_files = sorted(neutron_dir.glob('mesh_material_volumes*.npz'), + key=lambda p: int(p.stem.split('_')[-1]) + if p.stem[-1].isdigit() else 0) + if mmv_files: + self.results['mesh_material_volumes'] = [ + openmc.MeshMaterialVolumes.from_npz(f) for f in mmv_files + ] fluxes_file = neutron_dir / 'fluxes.npy' if fluxes_file.exists(): self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True)) @@ -665,10 +702,15 @@ def load_results(self, path: PathLike): # Load photon mesh material volumes if they exist (for mesh-based calculations) if self.method == 'mesh-based': - photon_mmv_file = photon_dir / 'mesh_material_volumes.npz' - if photon_mmv_file.exists(): - self.results['mesh_material_volumes_photon'] = \ - openmc.MeshMaterialVolumes.from_npz(photon_mmv_file) + photon_mmv_files = sorted( + photon_dir.glob('mesh_material_volumes*.npz'), + key=lambda p: int(p.stem.split('_')[-1]) + if p.stem[-1].isdigit() else 0) + if photon_mmv_files: + self.results['mesh_material_volumes_photon'] = [ + openmc.MeshMaterialVolumes.from_npz(f) + for f in photon_mmv_files + ] # Load tally IDs from JSON file tally_ids_path = photon_dir / 'tally_ids.json' diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 340ba616323..26529e6ce96 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -165,11 +165,11 @@ def capture_run(**kwargs): # Check that both tallies were created with the expected properties tally_names = [t.name for t in captured['tallies']] - assert 'MicroXS flux' in tally_names - assert 'MicroXS RR' in tally_names + assert 'MicroXS flux 0' in tally_names + assert 'MicroXS RR 0' in tally_names # Check that the RR tally has the expected nuclides and reactions - rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR') + rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR 0') assert rr.nuclides == ['U235'] assert rr.scores == ['fission'] diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py index a94f85c8c08..266bd49763a 100644 --- a/tests/unit_tests/test_r2s.py +++ b/tests/unit_tests/test_r2s.py @@ -6,7 +6,7 @@ @pytest.fixture -def simple_model_and_mesh(tmp_path): +def simple_model_and_mesh(): # Define two materials: water and Ni h2o = openmc.Material() h2o.add_nuclide("H1", 2.0) @@ -68,7 +68,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): nt = Path(outdir) / 'neutron_transport' assert (nt / 'fluxes.npy').exists() assert (nt / 'micros.h5').exists() - assert (nt / 'mesh_material_volumes.npz').exists() + assert (nt / 'mesh_material_volumes_0.npz').exists() act = Path(outdir) / 'activation' assert (act / 'depletion_results.h5').exists() pt = Path(outdir) / 'photon_transport' @@ -78,7 +78,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): # Basic results structure checks assert len(r2s.results['fluxes']) == 2 assert len(r2s.results['micros']) == 2 - assert len(r2s.results['mesh_material_volumes']) == 2 + assert len(r2s.results['mesh_material_volumes']) == 1 + assert len(r2s.results['mesh_material_volumes'][0]) == 2 assert len(r2s.results['activation_materials']) == 2 assert len(r2s.results['depletion_results']) == 2 @@ -93,11 +94,77 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): r2s_loaded.load_results(outdir) assert len(r2s_loaded.results['fluxes']) == 2 assert len(r2s_loaded.results['micros']) == 2 - assert len(r2s_loaded.results['mesh_material_volumes']) == 2 + assert len(r2s_loaded.results['mesh_material_volumes']) == 1 + assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2 assert len(r2s_loaded.results['activation_materials']) == 2 assert len(r2s_loaded.results['depletion_results']) == 2 +def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path): + model, _, _ = simple_model_and_mesh + + # Two 1x1x1 meshes that together cover the full domain, split along y. + # Each mesh element spans the full x range [-10, 10], crossing the x=0 + # material boundary, so both meshes contain both materials within their + # single element. + mesh1 = openmc.RegularMesh() + mesh1.lower_left = (-10.0, -10.0, -10.0) + mesh1.upper_right = (10.0, 0.0, 10.0) + mesh1.dimension = (1, 1, 1) + mesh2 = openmc.RegularMesh() + mesh2.lower_left = (-10.0, 0.0, -10.0) + mesh2.upper_right = (10.0, 10.0, 10.0) + mesh2.dimension = (1, 1, 1) + + r2s = R2SManager(model, [mesh1, mesh2]) + chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + + outdir = r2s.run( + timesteps=[(1.0, 'd')], + source_rates=[1.0], + photon_time_indices=[1], + output_dir=tmp_path, + chain_file=chain, + ) + + # Check that per-mesh MMV files were written + nt = Path(outdir) / 'neutron_transport' + assert (nt / 'fluxes.npy').exists() + assert (nt / 'micros.h5').exists() + assert (nt / 'mesh_material_volumes_0.npz').exists() + assert (nt / 'mesh_material_volumes_1.npz').exists() + act = Path(outdir) / 'activation' + assert (act / 'depletion_results.h5').exists() + pt = Path(outdir) / 'photon_transport' + assert (pt / 'tally_ids.json').exists() + assert (pt / 'time_1' / 'statepoint.10.h5').exists() + + # Two meshes, each with 1 element containing both materials → + # 2 element-material combinations per mesh, 4 total + assert len(r2s.results['mesh_material_volumes']) == 2 + assert len(r2s.results['mesh_material_volumes'][0]) == 2 + assert len(r2s.results['mesh_material_volumes'][1]) == 2 + assert len(r2s.results['fluxes']) == 4 + assert len(r2s.results['micros']) == 4 + assert len(r2s.results['activation_materials']) == 4 + assert len(r2s.results['depletion_results']) == 2 + + # Activation material names encode mesh index + amats = r2s.results['activation_materials'] + assert all(m.depletable for m in amats) + assert any('Mesh 0' in m.name for m in amats) + assert any('Mesh 1' in m.name for m in amats) + + # Check loading results + r2s_loaded = R2SManager(model, [mesh1, mesh2]) + r2s_loaded.load_results(outdir) + assert len(r2s_loaded.results['mesh_material_volumes']) == 2 + assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2 + assert len(r2s_loaded.results['mesh_material_volumes'][1]) == 2 + assert len(r2s_loaded.results['activation_materials']) == 4 + assert len(r2s_loaded.results['depletion_results']) == 2 + + def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): model, (c1, c2), _ = simple_model_and_mesh From c380716ef3861f48daf8e1419c009e0fa7459181 Mon Sep 17 00:00:00 2001 From: Jack Fletcher <115663563+j-fletcher@users.noreply.github.com> Date: Tue, 28 Apr 2026 17:10:03 -0400 Subject: [PATCH 07/11] Local adjoint source for Random Ray (#3717) --- docs/source/io_formats/settings.rst | 39 +- docs/source/methods/random_ray.rst | 48 +- docs/source/methods/variance_reduction.rst | 16 +- docs/source/usersguide/random_ray.rst | 55 +- docs/source/usersguide/variance_reduction.rst | 94 ++- .../openmc/random_ray/flat_source_domain.h | 8 +- .../openmc/random_ray/random_ray_simulation.h | 17 +- include/openmc/source.h | 1 + include/openmc/weight_windows.h | 3 + openmc/examples.py | 309 ++++++++ openmc/model/model.py | 44 ++ openmc/settings.py | 39 +- openmc/tallies.py | 91 ++- openmc/weight_windows.py | 57 +- src/random_ray/flat_source_domain.cpp | 101 ++- src/random_ray/random_ray_simulation.cpp | 217 +++--- src/settings.cpp | 22 +- src/source.cpp | 3 + src/weight_windows.cpp | 9 +- .../inputs_true.dat | 12 +- .../random_ray_adjoint_k_eff/inputs_true.dat | 12 +- .../random_ray_adjoint_local/__init__.py | 0 .../random_ray_adjoint_local/inputs_true.dat | 293 ++++++++ .../random_ray_adjoint_local/results_true.dat | 15 + .../random_ray_adjoint_local/test.py | 35 + .../infinite_medium/inputs_true.dat | 12 +- .../material_wise/inputs_true.dat | 12 +- .../stochastic_slab/inputs_true.dat | 12 +- .../infinite_medium/inputs_true.dat | 12 +- .../material_wise/inputs_true.dat | 12 +- .../stochastic_slab/inputs_true.dat | 12 +- .../infinite_medium/model/inputs_true.dat | 12 +- .../infinite_medium/user/inputs_true.dat | 12 +- .../stochastic_slab/model/inputs_true.dat | 12 +- .../stochastic_slab/user/inputs_true.dat | 12 +- .../infinite_medium/inputs_true.dat | 12 +- .../material_wise/inputs_true.dat | 12 +- .../stochastic_slab/inputs_true.dat | 12 +- .../eigen/inputs_true.dat | 12 +- .../fs/inputs_true.dat | 12 +- .../inputs_true.dat | 12 +- .../inputs_true.dat | 12 +- .../random_ray_entropy/settings.xml | 12 +- .../cell/inputs_true.dat | 12 +- .../material/inputs_true.dat | 12 +- .../universe/inputs_true.dat | 12 +- .../linear/inputs_true.dat | 12 +- .../linear_xy/inputs_true.dat | 12 +- .../flat/inputs_true.dat | 12 +- .../linear/inputs_true.dat | 12 +- .../False/inputs_true.dat | 12 +- .../True/inputs_true.dat | 12 +- .../flat/inputs_true.dat | 12 +- .../linear_xy/inputs_true.dat | 12 +- .../random_ray_halton_samples/inputs_true.dat | 12 +- .../random_ray_k_eff/inputs_true.dat | 12 +- .../random_ray_k_eff_mesh/inputs_true.dat | 12 +- .../random_ray_linear/linear/inputs_true.dat | 12 +- .../linear_xy/inputs_true.dat | 12 +- .../random_ray_low_density/inputs_true.dat | 12 +- .../inputs_true.dat | 12 +- .../random_ray_s2/inputs_true.dat | 12 +- .../random_ray_void/flat/inputs_true.dat | 12 +- .../random_ray_void/linear/inputs_true.dat | 12 +- .../hybrid/inputs_true.dat | 12 +- .../naive/inputs_true.dat | 12 +- .../simulation_averaged/inputs_true.dat | 12 +- .../hybrid/inputs_true.dat | 12 +- .../naive/inputs_true.dat | 12 +- .../simulation_averaged/inputs_true.dat | 12 +- .../weightwindows_fw_cadis/inputs_true.dat | 12 +- .../weightwindows_fw_cadis_local/__init__.py | 0 .../inputs_true.dat | 273 +++++++ .../results_true.dat | 696 ++++++++++++++++++ .../weightwindows_fw_cadis_local/test.py | 42 ++ .../flat/inputs_true.dat | 12 +- .../linear/inputs_true.dat | 12 +- 77 files changed, 2664 insertions(+), 463 deletions(-) create mode 100644 tests/regression_tests/random_ray_adjoint_local/__init__.py create mode 100644 tests/regression_tests/random_ray_adjoint_local/inputs_true.dat create mode 100644 tests/regression_tests/random_ray_adjoint_local/results_true.dat create mode 100644 tests/regression_tests/random_ray_adjoint_local/test.py create mode 100644 tests/regression_tests/weightwindows_fw_cadis_local/__init__.py create mode 100644 tests/regression_tests/weightwindows_fw_cadis_local/inputs_true.dat create mode 100644 tests/regression_tests/weightwindows_fw_cadis_local/results_true.dat create mode 100644 tests/regression_tests/weightwindows_fw_cadis_local/test.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index a50922b041e..d63376e2bdc 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -597,7 +597,7 @@ found in the :ref:`random ray user guide `. *Default*: None - :source: + :ray_source: Specifies the starting ray distribution, and follows the format for :ref:`source_element`. It must be uniform in space and angle and cover the full domain. It does not represent a physical neutron or photon source -- it @@ -605,6 +605,35 @@ found in the :ref:`random ray user guide `. *Default*: None + :adjoint_source: + Specifies an adjoint fixed source for adjoint transport simulations, and + follows the format for :ref:`source_element`. The distributions which make + up the adjoint source are subject to the same restrictions as forward + fixed sources in Random Ray mode. + + *Default*: None + + :adjoint: + Specifies whether to perform adjoint transport. The default is 'False', + corresponding to forward transport. + + *Default*: None + + :volume_estimator: + Specifies choice of volume estimator for the random ray solver. Options + are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'. + + *Default*: None + + :volume_normalized_flux_tallies: + Specifies whether to normalize flux tallies by volume (bool). The + default is 'False'. When enabled, flux tallies will be reported in units + of cm/cm^3. When disabled, flux tallies will be reported in units of cm + (i.e., total distance traveled by neutrons in the spatial tally + region). + + *Default*: None + :sample_method: Specifies the method for sampling the starting ray distribution. This element can be set to "prng" or "halton". @@ -1696,6 +1725,14 @@ mesh-based weight windows. The ratio of the lower to upper weight window bounds. *Default*: 5.0 + + For FW-CADIS: + + :targets: + A sequence of IDs corresponding to the tallies which cover phase + space regions of interest for local variance reduction. + + *Default*: None --------------------------------------- ```` Element diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 5e17316aa1c..8bc2a0a1bf5 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -1081,28 +1081,32 @@ lifetimes. In OpenMC, the random ray adjoint solver is implemented simply by transposing the scattering matrix, swapping :math:`\nu\Sigma_f` and :math:`\chi`, and then -running a normal transport solve. When no external fixed source is present, no -additional changes are needed in the transport process. However, if an external -fixed forward source is present in the simulation problem, then an additional -step is taken to compute the accompanying fixed adjoint source. In OpenMC, the -adjoint flux does *not* represent a response function for a particular detector -region. Rather, the adjoint flux is the global response, making it appropriate -for use with weight window generation schemes for global variance reduction. -Thus, if using a fixed source, the external source for the adjoint mode is -simply computed as being :math:`1 / \phi`, where :math:`\phi` is the forward -scalar flux that results from a normal forward solve (which OpenMC will run -first automatically when in adjoint mode). The adjoint external source will be -computed for each source region in the simulation mesh, independent of any -tallies. The adjoint external source is always flat, even when a linear -scattering and fission source shape is used. When in adjoint mode, all reported -results (e.g., tallies, eigenvalues, etc.) are derived from the adjoint flux, -even when the physical meaning is not necessarily obvious. These values are -still reported, though we emphasize that the primary use case for adjoint mode -is for producing adjoint flux tallies to support subsequent perturbation studies -and weight window generation. - -Note that the adjoint :math:`k_{eff}` is statistically the same as the forward -:math:`k_{eff}`, despite the flux distributions taking different shapes. +running a normal transport solve. When no external fixed forward source is +present, or if an adjoint fixed source is specifically provided, no additional +changes are needed in the transport process. This adjoint source can +correspond, for example, to a detector response function in a particular +region. However, if an external fixed forward source is present in the +simulation problem without an adjoint fixed source, an additional step is taken +to compute the accompanying forward-weighted adjoint source. In this case, the +adjoint flux does *not* represent the importance of locations in phase space to +detector response; rather, the "response" in question is a uniform distribution +of Monte Carlo particle density, making the importance provided by the adjoint +flux appropriate for use with weight window generation schemes for global +variance reduction. Thus, if using a fixed source, the forward-weighted +external source for adjoint mode is simply computed as being :math:`1 / \phi`, +where :math:`\phi` is the forward scalar flux that results from a normal +forward solve (which OpenMC will run first automatically when in adjoint mode). +The adjoint external source will be computed for each source region in the +simulation mesh, independent of any tallies. The adjoint external source is +always flat, even when a linear scattering and fission source shape is used. + +When in adjoint mode, all reported results (e.g., tallies, eigenvalues, etc.) +are derived from the adjoint flux, even when the physical meaning is not +necessarily obvious. These values are still reported, though we emphasize that +the primary use case for adjoint mode is for producing adjoint flux tallies to +support subsequent perturbation studies and weight window generation. Note +however that the adjoint :math:`k_{eff}` is statistically the same as the +forward :math:`k_{eff}`, despite the flux distributions taking different shapes. --------------------------- Fundamental Sources of Bias diff --git a/docs/source/methods/variance_reduction.rst b/docs/source/methods/variance_reduction.rst index cdda5ea9295..7778e0714ce 100644 --- a/docs/source/methods/variance_reduction.rst +++ b/docs/source/methods/variance_reduction.rst @@ -82,8 +82,8 @@ where it was born from. The Forward-Weighted Consistent Adjoint Driven Importance Sampling method, or `FW-CADIS method `_, produces weight windows -for global variance reduction given adjoint flux information throughout the -entire domain. The weight window lower bound is defined in Equation +for global or local variance reduction given adjoint flux information throughout +the entire domain. The weight window lower bound is defined in Equation :eq:`fw_cadis`, and also involves a normalization step not shown here. .. math:: @@ -135,6 +135,18 @@ aware of this. \text{FOM} = \frac{1}{\text{Time} \times \sigma^2} +Finally, one unique capability of the FW-CADIS weight window generator is to +produce weight windows for local variance reduction, given a list of the +responses of interest. This is controlled by optionally specifying target +tallies from the :class:`openmc.model.Model` to the +:class:`openmc.WeightWindowGenerator`, as illustrated in the +:ref:`user guide`. If target tallies for local variance +reduction are supplied, then the adjoint sources are only populated after the +initial forward simulation in the source regions associated with those tallies. +In other regions, the adjoint source term is instead set to zero. The Random +Ray solver then determines the adjoint flux map used to generate FW-CADIS +weight windows following the usual technique. + .. _methods_source_biasing: -------------- diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 0c9a0402814..d35aff83b78 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -944,6 +944,8 @@ as:: which will greatly improve the quality of the linear source term in 2D simulations. +.. _usersguide_random_ray_run_modes: + --------------------------------- Fixed Source and Eigenvalue Modes --------------------------------- @@ -1073,22 +1075,47 @@ The adjoint flux random ray solver mode can be enabled as:: settings.random_ray['adjoint'] = True -When enabled, OpenMC will first run a forward transport simulation followed by -an adjoint transport simulation. The purpose of the forward solve is to compute -the adjoint external source when an external source is present in the -simulation. Simulation settings (e.g., number of rays, batches, etc.) will be -identical for both simulations. At the conclusion of the run, all results (e.g., -tallies, plots, etc.) will be derived from the adjoint flux rather than the -forward flux but are not labeled any differently. The initial forward flux -solution will not be stored or available in the final statepoint file. Those -wishing to do analysis requiring both the forward and adjoint solutions will -need to run two separate simulations and load both statepoint files. +When enabled, OpenMC will first run a forward transport simulation if there are +no user-specified adjoint sources present, followed by an adjoint transport +simulation. Fixed adjoint sources can be specified on the +:attr:`openmc.Settings.random_ray` dictionary as follows:: + + # Geometry definition + ... + detector_cell = openmc.Cell(fill=detector_mat, name='cell where detector will be') + ... + # Define fixed adjoint neutron source + strengths = [1.0] + midpoints = [1.0e-4] + energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths) + + adj_source = openmc.IndependentSource( + energy=energy_distribution, + constraints={'domains': [detector_cell]} + ) + + # Add to random_ray dict + settings.random_ray['adjoint_source'] = adj_source + +The same constraints apply to the user-defined adjoint source as to the forward +source, described in the :ref:`Fixed Source and Eigenvalue section +`. If this source is not provided, a forward +solve must take place to compute the adjoint external source when a forward +external source is present in the problem. Simulation settings (e.g., number of +rays, batches, etc.) will be identical for both calculations. At the +conclusion of the run, all results (e.g., tallies, plots, etc.) will be +derived from the adjoint flux rather than the forward flux but are not labeled +any differently. The initial forward flux solution will not be stored or +available in the final statepoint file. Those wishing to do analysis requiring +both the forward and adjoint solutions will need to run two separate +simulations and load both statepoint files. .. note:: - When adjoint mode is selected, OpenMC will always perform a full forward - solve and then run a full adjoint solve immediately afterwards. Statepoint - and tally results will be derived from the adjoint flux, but will not be - labeled any differently. + Use of the automated + :ref:`FW-CADIS weight window generator` is not + currently compatible with user-defined adjoint sources. Instead, the + initial forward calculation is used to assign "forward-weighted" adjoint + sources to the tally regions of interest. --------------------------------------- Putting it All Together: Example Inputs diff --git a/docs/source/usersguide/variance_reduction.rst b/docs/source/usersguide/variance_reduction.rst index 8d41807e1de..d551195f5cc 100644 --- a/docs/source/usersguide/variance_reduction.rst +++ b/docs/source/usersguide/variance_reduction.rst @@ -4,26 +4,27 @@ Variance Reduction ================== -Global variance reduction in OpenMC is accomplished by weight windowing -or source biasing techniques, the latter of which additionally provides a -local variance reduction capability. OpenMC is capable of generating weight -windows using either the MAGIC or FW-CADIS methods. Both techniques will -produce a ``weight_windows.h5`` file that can be loaded and used later on. In +Global and local variance reduction are possible in OpenMC through both weight +windowing and source biasing techniques. OpenMC is capable of generating weight +windows using either the MAGIC or FW-CADIS methods, the latter with an optional +capability for local variance reduction. Both techniques will produce a +``weight_windows.h5`` file that can be loaded and used later on. In this section, we first break down the steps required to generate and apply weight windows, then describe how source biasing may be applied. .. _ww_generator: ------------------------------------- -Generating Weight Windows with MAGIC ------------------------------------- +------------------------------------------- +Generating Global Weight Windows with MAGIC +------------------------------------------- As discussed in the :ref:`methods section `, MAGIC is an iterative method that uses flux tally information from a Monte Carlo -simulation to produce weight windows for a user-defined mesh. While generating -the weight windows, OpenMC is capable of applying the weight windows generated -from a previous batch while processing the next batch, allowing for progressive -improvement in the weight window quality across iterations. +simulation to produce weight windows for a user-defined mesh with the objective +of global variance reduction. While generating the weight windows, OpenMC is +capable of applying the weight windows generated from a previous batch while +processing the next batch, allowing for progressive improvement in the weight +window quality across iterations. The typical way of generating weight windows is to define a mesh and then add an :class:`openmc.WeightWindowGenerator` object to an :attr:`openmc.Settings` @@ -71,15 +72,20 @@ At the end of the simulation, a ``weight_windows.h5`` file will be saved to disk for later use. Loading it in another subsequent simulation will be discussed in the "Using Weight Windows" section below. ------------------------------------------------------- -Generating Weight Windows with FW-CADIS and Random Ray ------------------------------------------------------- +.. _usersguide_fw_cadis: + +---------------------------------------------------------------------- +Generating Global or Local Weight Windows with FW-CADIS and Random Ray +---------------------------------------------------------------------- Weight window generation with FW-CADIS and random ray in OpenMC uses the same -exact strategy as with MAGIC. An :class:`openmc.WeightWindowGenerator` object is -added to the :attr:`openmc.Settings` object, and a ``weight_windows.h5`` will be -generated at the end of the simulation. The only difference is that the code -must be run in random ray mode. A full description of how to enable and setup +exact strategy as with MAGIC. Using FW-CADIS, however, also enables +local variance reduction in fixed source problems through the :attr:`targets` +attribute, which is described later in this section. To enable FW-CADIS, an +:class:`openmc.WeightWindowGenerator` object is added to the +:attr:`openmc.Settings` object, and a ``weight_windows.h5`` will be generated +at the end of the simulation. The only procedural difference is that the code +must be run in random ray mode. A full description of how to enable and setup random ray mode can be found in the :ref:`Random Ray User Guide `. .. note:: @@ -90,7 +96,7 @@ random ray mode can be found in the :ref:`Random Ray User Guide `. ray solver. A high level overview of the current workflow for generation of weight windows with FW-CADIS using random ray is given below. -1. Begin by making a deepy copy of your continuous energy Python model and then +1. Begin by making a deep copy of your continuous energy Python model and then convert the copy to be multigroup and use the random ray transport solver. The conversion process can largely be automated as described in more detail in the :ref:`random ray quick start guide `, summarized below:: @@ -148,7 +154,53 @@ random ray mode can be found in the :ref:`Random Ray User Guide `. assigning to ``model.settings.random_ray['source_region_meshes']``) and for weight window generation. -3. When running your multigroup random ray input deck, OpenMC will automatically +3. (Optional) If local variance reduction is desired in a fixed-source problem, + populate the :attr:`targets` attribute with an :class:`openmc.Tallies` + instance or an iterable of tally IDs indicating the tallies of interest for + variance reduction:: + + # Build a new example and WWG for local variance reduction + from openmc.examples import random_ray_three_region_cube_with_detectors + new_model = random_ray_three_region_cube_with_detectors() + + ww_mesh = openmc.RegularMesh() + n = 7 + width = 35.0 + ww_mesh.dimension = (n, n, n) + ww_mesh.lower_left = (0.0, 0.0, 0.0) + ww_mesh.upper_right = (width, width, width) + + wwg = openmc.WeightWindowGenerator( + method="fw_cadis", + mesh=ww_mesh, + max_realizations=new_model.settings.batches + ) + new_model.settings.weight_window_generators = wwg + new_model.settings.random_ray['volume_estimator'] = 'naive' + + # Get the tallies of interest + target_tallies = openmc.Tallies() + + for tally in list(new_model.tallies): + if tally.name in {"Detector 1 Tally", "Detector 2 Tally"}: + target_tallies.append(tally) + + # Add to WeightWindowGenerator + wwg.targets = target_tallies + +.. warning:: + The tallies designated as FW-CADIS targets to the + :class:`~openmc.WeightWindowGenerator` must be present under the + :class:`~openmc.model.Model.tallies` attribute of the + :class:`~openmc.model.Model` as well in order to be recognized as valid + local variance reduction targets. This check is performed when the + :func:`openmc.model.Model.export_to_model_xml` or + :func:`openmc.model.Model.export_to_xml` functions are called, meaning + that the standalone :func:`openmc.Settings.export_to_xml` and + :func:`openmc.Tallies.export_to_xml` methods should not be used with + FW-CADIS local variance reduction. + +4. When running your multigroup random ray input deck, OpenMC will automatically run a forward solve followed by an adjoint solve, with a ``weight_windows.h5`` file generated at the end. The ``weight_windows.h5`` file will contain FW-CADIS generated weight windows. This file can be used in diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index c40982712eb..6f51af34ded 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -40,9 +40,10 @@ class FlatSourceDomain { void random_ray_tally(); virtual void accumulate_iteration_flux(); void output_to_vtk() const; - void convert_external_sources(); + void convert_external_sources(bool use_adjoint_sources); void count_external_source_regions(); - void set_adjoint_sources(); + void set_fw_adjoint_sources(); + void set_local_adjoint_sources(); void flux_swap(); virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; @@ -76,6 +77,7 @@ class FlatSourceDomain { // Static Data members static bool volume_normalized_flux_tallies_; static bool adjoint_; // If the user wants outputs based on the adjoint flux + static bool fw_cadis_local_; static double diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization // for transport corrected MGXS data @@ -84,6 +86,8 @@ class FlatSourceDomain { static std::unordered_map>> mesh_domain_map_; + static std::vector fw_cadis_local_targets_; + //---------------------------------------------------------------------------- // Static data members static RandomRayVolumeEstimator volume_estimator_; diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index 68c7779edad..ccd2cbe4763 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -20,8 +20,9 @@ class RandomRaySimulation { //---------------------------------------------------------------------------- // Methods void apply_fixed_sources_and_mesh_domains(); - void prepare_fixed_sources_adjoint(); - void prepare_adjoint_simulation(); + void prepare_fw_fixed_sources_adjoint(); + void prepare_local_fixed_sources_adjoint(); + void prepare_adjoint_simulation(bool fw_adjoint); void simulate(); void output_simulation_results() const; void instability_check( @@ -34,15 +35,9 @@ class RandomRaySimulation { // Accessors FlatSourceDomain* domain() const { return domain_.get(); } - //---------------------------------------------------------------------------- - // Public data members - - // Flag for adjoint simulation; - bool adjoint_needed_; - private: //---------------------------------------------------------------------------- - // Private data members + // Data members // Contains all flat source region data unique_ptr domain_; @@ -57,9 +52,6 @@ class RandomRaySimulation { // Number of energy groups int negroups_; - // Toggle for first simulation - bool is_first_simulation_; - }; // class RandomRaySimulation //============================================================================ @@ -67,7 +59,6 @@ class RandomRaySimulation { //============================================================================ void validate_random_ray_inputs(); -void print_adjoint_header(); void openmc_finalize_random_ray(); } // namespace openmc diff --git a/include/openmc/source.h b/include/openmc/source.h index 1ef7eba2b24..e307b1ed217 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -43,6 +43,7 @@ class Source; namespace model { extern vector> external_sources; +extern vector> adjoint_sources; // Probability distribution for selecting external sources extern DiscreteIndex external_sources_probability; diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 42846f9d140..a5d404133ce 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -223,6 +223,9 @@ class WeightWindowsGenerator { double threshold_ {1.0}; // targets_; }; //============================================================================== diff --git a/openmc/examples.py b/openmc/examples.py index f7f2bd48da1..90a0bffe7b8 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -1310,3 +1310,312 @@ def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): model.tallies = tallies return model + +def random_ray_three_region_cube_with_detectors() -> openmc.Model: + """Create a three region cube model with two external tally regions. + + This is an adaptation of the simple monoenergetic problem of a cube with + three concentric cubic regions. The innermost region is near void (with + Sigma_t around 10^-5) and contains an external isotropic source term, the + middle region is a mild scatterer (with Sigma_t around 10^-3), and the + outer region of the cube is a scatterer and absorber (with Sigma_t around + 1). + + Two cubic "detector" regions are found outside this geometry, one along the + y-axis near z=0, and the other in the upper right corner of the system. + The size of each detector is scaled to be equal to that of the source + region. The model returned by this function contains cell tallies on each + detector. + + Returns + ------- + model : openmc.Model + A three region cube model + + """ + + model = openmc.Model() + + ########################################################################### + # Helper function creates a 3 region cube with different fills in each region + def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): + cube = [[[0 for _ in range(N)] for _ in range(N)] for _ in range(N)] + for i in range(N): + for j in range(N): + for k in range(N): + if i < n_1 and j >= (N-n_1) and k < n_1: + cube[i][j][k] = fill_1 + elif i < n_2 and j >= (N-n_2) and k < n_2: + cube[i][j][k] = fill_2 + else: + cube[i][j][k] = fill_3 + return cube + + ########################################################################### + # Create multigroup data + + # Instantiate the energy group data + ebins = [1e-5, 20.0e6] + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + cavity_sigma_a = 4.0e-5 + cavity_sigma_s = 3.0e-3 + cavity_mat_data = openmc.XSdata('cavity', groups) + cavity_mat_data.order = 0 + cavity_mat_data.set_total([cavity_sigma_a + cavity_sigma_s]) + cavity_mat_data.set_absorption([cavity_sigma_a]) + cavity_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[cavity_sigma_s]]]), 0, 3)) + + absorber_sigma_a = 0.50 + absorber_sigma_s = 0.50 + absorber_mat_data = openmc.XSdata('absorber', groups) + absorber_mat_data.order = 0 + absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s]) + absorber_mat_data.set_absorption([absorber_sigma_a]) + absorber_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3)) + + multiplier = 0.01 + source_sigma_a = cavity_sigma_a * multiplier + source_sigma_s = cavity_sigma_s * multiplier + source_mat_data = openmc.XSdata('source', groups) + source_mat_data.order = 0 + source_mat_data.set_total([source_sigma_a + source_sigma_s]) + source_mat_data.set_absorption([source_sigma_a]) + source_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3)) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas( + [source_mat_data, cavity_mat_data, absorber_mat_data]) + mg_cross_sections_file.export_to_hdf5() + + ########################################################################### + # Create materials for the problem + + # Instantiate some Macroscopic Data + source_data = openmc.Macroscopic('source') + cavity_data = openmc.Macroscopic('cavity') + absorber_data = openmc.Macroscopic('absorber') + + # Instantiate some Materials and register the appropriate Macroscopic objects + source_mat = openmc.Material(name='source') + source_mat.set_density('macro', 1.0) + source_mat.add_macroscopic(source_data) + + cavity_mat = openmc.Material(name='cavity') + cavity_mat.set_density('macro', 1.0) + cavity_mat.add_macroscopic(cavity_data) + + absorber_mat = openmc.Material(name='absorber') + absorber_mat.set_density('macro', 1.0) + absorber_mat.add_macroscopic(absorber_data) + + # Instantiate a Materials collection + materials_file = openmc.Materials([source_mat, cavity_mat, absorber_mat]) + materials_file.cross_sections = "mgxs.h5" + + ########################################################################### + # Define problem geometry + + source_cell = openmc.Cell(fill=source_mat, name='infinite source region') + cavity_cell = openmc.Cell(fill=cavity_mat, name='cube cavity region') + absorber_cell = openmc.Cell( + fill=absorber_mat, name='absorber region') + + source_universe = openmc.Universe(name='source universe') + source_universe.add_cells([source_cell]) + + cavity_universe = openmc.Universe() + cavity_universe.add_cells([cavity_cell]) + + absorber_universe = openmc.Universe() + absorber_universe.add_cells([absorber_cell]) + + absorber_width = 30.0 + n_base = 6 + + # This variable can be increased above 1 to refine the FSR mesh resolution further + refinement_level = 2 + + n = n_base * refinement_level + pitch = absorber_width / n + + pattern = fill_cube(n, 1*refinement_level, 5*refinement_level, + source_universe, cavity_universe, absorber_universe) + + lattice = openmc.RectLattice() + lattice.lower_left = [0.0, 0.0, 0.0] + lattice.pitch = [pitch, pitch, pitch] + lattice.universes = pattern + + lattice_cell = openmc.Cell(fill=lattice) + + lattice_uni = openmc.Universe() + lattice_uni.add_cells([lattice_cell]) + + x_low = openmc.XPlane(x0=0.0, boundary_type='reflective') + x_high = openmc.XPlane(x0=absorber_width) + + y_low = openmc.YPlane(y0=0.0, boundary_type='reflective') + y_high = openmc.YPlane(y0=absorber_width) + + z_low = openmc.ZPlane(z0=0.0, boundary_type='reflective') + z_high = openmc.ZPlane(z0=absorber_width) + + cube_domain = openmc.Cell(fill=lattice_uni, region=+x_low & - + x_high & +y_low & -y_high & +z_low & -z_high, name='full domain') + + detect_width = absorber_width / n_base + outer_width = absorber_width + detect_width + + x_outer = openmc.XPlane(x0=outer_width, boundary_type='vacuum') + y_outer = openmc.YPlane(y0=outer_width, boundary_type='vacuum') + z_outer = openmc.ZPlane(z0=outer_width, boundary_type='vacuum') + + detector1_right = openmc.XPlane(x0=detect_width) + detector1_top = openmc.ZPlane(z0=detect_width) + + detector1_region = ( + +x_low & -detector1_right & + +y_high & -y_outer & + +z_low & -detector1_top + ) + detector1 = openmc.Cell( + name='detector 1', + fill=absorber_mat, + region=detector1_region + ) + + detector2_region = ( + +x_high & -x_outer & + +y_high & -y_outer & + +z_high & -z_outer + ) + detector2 = openmc.Cell( + name='detector 2', + fill=absorber_mat, + region=detector2_region + ) + + external_x = ( + +x_high & +y_low & +z_low & -x_outer & + ((-y_outer & -z_high) | (-y_high & +z_high & -z_outer)) + ) + external_y = ( + +y_high & -y_outer & + ( + (+detector1_right & -x_high & +z_low & -z_outer) | + (-detector1_right & +x_low & +detector1_top & -z_outer) | + (+x_high & -x_outer & +z_low & -z_high) + ) + ) + external_z = ( + +x_low & +y_low & +z_high & -z_outer & + ((-y_outer & -x_high) | (-y_high & +x_high & -x_outer)) + ) + external_cell = openmc.Cell(fill=cavity_mat, + region=(external_x | external_y | external_z), + name='outside cube') + + root = openmc.Universe( + name='root universe', + cells=[cube_domain, detector1, detector2, external_cell] + ) + + # Create a geometry with the two cells and export to XML + geometry = openmc.Geometry(root) + + ########################################################################### + # Define problem settings + + # Instantiate a Settings object, set all runtime parameters, and export to XML + settings = openmc.Settings() + settings.energy_mode = "multi-group" + settings.inactive = 5 + settings.batches = 10 + settings.particles = 500 + settings.run_mode = 'fixed source' + + # Create an initial uniform spatial source for ray integration + lower_left_ray = [0.0, 0.0, 0.0] + upper_right_ray = [outer_width, outer_width, outer_width] + uniform_dist_ray = openmc.stats.Box( + lower_left_ray, upper_right_ray, only_fissionable=False) + rr_source = openmc.IndependentSource(space=uniform_dist_ray) + + settings.random_ray['distance_active'] = 800.0 + settings.random_ray['distance_inactive'] = 100.0 + settings.random_ray['ray_source'] = rr_source + settings.random_ray['volume_normalized_flux_tallies'] = True + + # Create a rectilinear source region mesh + sr_mesh = openmc.RegularMesh() + sr_mesh.dimension = (14, 14, 14) + sr_mesh.lower_left = (0.0, 0.0, 0.0) + sr_mesh.upper_right = (outer_width, outer_width, outer_width) + settings.random_ray['source_region_meshes'] = [(sr_mesh, [root])] + + # Create the neutron source in the bottom right of the moderator + # Good - fast group appears largest (besides most thermal) + strengths = [1.0] + midpoints = [100.0] + energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths) + + source = openmc.IndependentSource(energy=energy_distribution, constraints={ + 'domains': [source_universe]}, strength=3.14) + + settings.source = [source] + + ########################################################################### + # Define tallies + + estimator = 'tracklength' + + detector1_filter = openmc.CellFilter(detector1) + detector1_tally = openmc.Tally(name="Detector 1 Tally") + detector1_tally.filters = [detector1_filter] + detector1_tally.scores = ['flux'] + detector1_tally.estimator = estimator + + detector2_filter = openmc.CellFilter(detector2) + detector2_tally = openmc.Tally(name="Detector 2 Tally") + detector2_tally.filters = [detector2_filter] + detector2_tally.scores = ['flux'] + detector2_tally.estimator = estimator + + absorber_filter = openmc.MaterialFilter(absorber_mat) + absorber_tally = openmc.Tally(name="Absorber Tally") + absorber_tally.filters = [absorber_filter] + absorber_tally.scores = ['flux'] + absorber_tally.estimator = estimator + + cavity_filter = openmc.MaterialFilter(cavity_mat) + cavity_tally = openmc.Tally(name="Cavity Tally") + cavity_tally.filters = [cavity_filter] + cavity_tally.scores = ['flux'] + cavity_tally.estimator = estimator + + source_filter = openmc.MaterialFilter(source_mat) + source_tally = openmc.Tally(name="Source Tally") + source_tally.filters = [source_filter] + source_tally.scores = ['flux'] + source_tally.estimator = estimator + + # Instantiate a Tallies collection and export to XML + tallies = openmc.Tallies([detector1_tally, + detector2_tally, + absorber_tally, + cavity_tally, + source_tally]) + + ########################################################################### + # Assmble Model + + model.geometry = geometry + model.materials = materials_file + model.settings = settings + model.tallies = tallies + + return model \ No newline at end of file diff --git a/openmc/model/model.py b/openmc/model/model.py index 884e3ff6f95..c9a24b8b389 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -265,6 +265,46 @@ def add_kinetics_parameters_tallies(self, num_groups: int | None = None): denom_tally = openmc.Tally(name='IFP denominator') denom_tally.scores = ['ifp-denominator'] self.tallies.append(denom_tally) + + # TODO: This should also be incorporated into lower-level calls in + # settings.py, but it requires information about the tallies currently + # on the active Model + def _assign_fw_cadis_tally_IDs(self): + # Verify that all tallies assigned as targets on WeightWindowGenerators + # exist within model.tallies. If this is the case, convert the .targets + # attribute of each WeightWindowGenerator to a sequence of tally IDs. + if len(self.settings.weight_window_generators) == 0: + return + + # List of valid tally IDs + reference_tally_ids = np.asarray([tal.id for tal in self.tallies]) + + for wwg in self.settings.weight_window_generators: + # Only proceeds if the "targets" attribute is an openmc.Tallies, + # which means it hasn't been checked against model.tallies. + if isinstance(wwg.targets, openmc.Tallies): + id_vec = [] + for tal in wwg.targets: + # check against model tallies for equivalence + id_next = None + for reference_tal in self.tallies: + if tal == reference_tal: + id_next = reference_tal.id + break + + if id_next == None: + raise RuntimeError( + f'Local FW-CADIS target tally {tal.id} not found on model.tallies!') + else: + id_vec.append(id_next) + + wwg.targets = id_vec + + elif isinstance(wwg.targets, np.ndarray): + invalid = wwg.targets[~np.isin(wwg.targets, reference_tally_ids)] + if len(invalid) > 0: + raise RuntimeError( + f'Local FW-CADIS target tally IDs {invalid} not found on model.tallies!') @classmethod def from_xml( @@ -576,6 +616,7 @@ def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False, if not d.is_dir(): d.mkdir(parents=True, exist_ok=True) + self._assign_fw_cadis_tally_IDs() self.settings.export_to_xml(d) self.geometry.export_to_xml(d, remove_surfs=remove_surfs) @@ -634,6 +675,9 @@ def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool = "set the Geometry.merge_surfaces attribute instead.") self.geometry.merge_surfaces = True + # Link FW-CADIS WeightWindowGenerator target tallies, if present + self._assign_fw_cadis_tally_IDs() + # provide a memo to track which meshes have been written mesh_memo = set() settings_element = self.settings.to_xml_element(mesh_memo) diff --git a/openmc/settings.py b/openmc/settings.py index 6919afca4cc..1090babda49 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -236,6 +236,9 @@ class Settings: stabilization, which may be desirable as stronger diagonal stabilization also tends to dampen the convergence rate of the solver, thus requiring more iterations to converge. + :adjoint_source: + Source object used to define localized adjoint source/detector response + function. .. versionadded:: 0.15.0 resonance_scattering : dict @@ -1421,6 +1424,14 @@ def random_ray(self, random_ray: dict): cv.check_type('diagonal stabilization rho', value, Real) cv.check_greater_than('diagonal stabilization rho', value, 0.0, True) + elif key == 'adjoint_source': + if not isinstance(value, MutableSequence): + value = [value] + for source in value: + if not isinstance(source, SourceBase): + raise ValueError( + f'Invalid adjoint source type: {type(source)}. ' + 'Expected openmc.SourceBase.') else: raise ValueError(f'Unable to set random ray to "{key}" which is ' 'unsupported by OpenMC') @@ -1973,11 +1984,12 @@ def _create_random_ray_subelement(self, root, mesh_memo=None): element = ET.SubElement(root, "random_ray") for key, value in self._random_ray.items(): if key == 'ray_source' and isinstance(value, SourceBase): + subelement = ET.SubElement(element, 'ray_source') source_element = value.to_xml_element() if source_element.find('bias') is not None: raise RuntimeError( "Ray source distributions should not be biased.") - element.append(source_element) + subelement.append(source_element) elif key == 'source_region_meshes': subelement = ET.SubElement(element, 'source_region_meshes') @@ -1995,8 +2007,20 @@ def _create_random_ray_subelement(self, root, mesh_memo=None): path = f"./mesh[@id='{mesh.id}']" if root.find(path) is None: root.append(mesh.to_xml_element()) - if mesh_memo is not None: + if mesh_memo is not None: mesh_memo.add(mesh.id) + elif key == 'adjoint_source': + subelement = ET.SubElement(element, 'adjoint_source') + # Check that all entries are valid SourceBase instances, in case + # the random_ray setter was not used to populate dict entries. + if not isinstance(value, MutableSequence): + value = [value] + for source in value: + if not isinstance(source, SourceBase): + raise ValueError( + f'Invalid adjoint source type: {type(source)}. ' + 'Expected openmc.SourceBase.') + subelement.append(source.to_xml_element()) elif isinstance(value, bool): subelement = ET.SubElement(element, key) subelement.text = str(value).lower() @@ -2443,8 +2467,9 @@ def _random_ray_from_xml_element(self, root, meshes=None): for child in elem: if child.tag in ('distance_inactive', 'distance_active', 'diagonal_stabilization_rho'): self.random_ray[child.tag] = float(child.text) - elif child.tag == 'source': - source = SourceBase.from_xml_element(child) + elif child.tag == 'ray_source': + source_element = child.find('source') + source = SourceBase.from_xml_element(source_element) if child.find('bias') is not None: raise RuntimeError( "Ray source distributions should not be biased.") @@ -2461,6 +2486,12 @@ def _random_ray_from_xml_element(self, root, meshes=None): self.random_ray['adjoint'] = ( child.text in ('true', '1') ) + elif child.tag == 'adjoint_source': + self.random_ray['adjoint_source'] = [] + for subelem in child.findall('source'): + src = SourceBase.from_xml_element(subelem) + # add newly constructed source object to the list + self.random_ray['adjoint_source'].append(src) elif child.tag == 'sample_method': self.random_ray['sample_method'] = child.text elif child.tag == 'source_region_meshes': diff --git a/openmc/tallies.py b/openmc/tallies.py index 151add2bee3..ceced425533 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -16,6 +16,23 @@ import openmc import openmc.checkvalue as cv +from openmc.filter import ( + Filter, + DistribcellFilter, + EnergyFunctionFilter, + DelayedGroupFilter, + FilterMeta, + MeshFilter, + MeshBornFilter, +) +from openmc.arithmetic import ( + CrossFilter, + AggregateFilter, + CrossScore, + AggregateScore, + CrossNuclide, + AggregateNuclide, +) from ._sparse_compat import lil_array from ._xml import clean_indentation, get_elem_list, get_text from .mixin import IDManagerMixin @@ -31,9 +48,9 @@ # The following indicate acceptable types when setting Tally.scores, # Tally.nuclides, and Tally.filters -_SCORE_CLASSES = (str, openmc.CrossScore, openmc.AggregateScore) -_NUCLIDE_CLASSES = (str, openmc.CrossNuclide, openmc.AggregateNuclide) -_FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter) +_SCORE_CLASSES = (str, CrossScore, AggregateScore) +_NUCLIDE_CLASSES = (str, CrossNuclide, AggregateNuclide) +_FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter) # Valid types of estimators ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'} @@ -421,7 +438,7 @@ def _read_results(self): self._num_realizations = int(group['n_realizations'][()]) for filt in self.filters: - if isinstance(filt, openmc.DistribcellFilter): + if isinstance(filt, DistribcellFilter): filter_group = f[f'tallies/filters/filter {filt.id}'] filt._num_bins = int(filter_group['n_bins'][()]) @@ -1089,8 +1106,8 @@ def _can_merge_filters(self, other): return False # Return False if only one tally has a delayed group filter - tally1_dg = self.contains_filter(openmc.DelayedGroupFilter) - tally2_dg = other.contains_filter(openmc.DelayedGroupFilter) + tally1_dg = self.contains_filter(DelayedGroupFilter) + tally2_dg = other.contains_filter(DelayedGroupFilter) if tally1_dg != tally2_dg: return False @@ -1602,7 +1619,7 @@ def find_filter(self, filter_type): # Also check to see if the desired filter is wrapped up in an # aggregate - elif isinstance(test_filter, openmc.AggregateFilter): + elif isinstance(test_filter, AggregateFilter): if isinstance(test_filter.aggregate_filter, filter_type): return test_filter @@ -1704,7 +1721,7 @@ def get_filter_indices(self, filters=[], filter_bins=[]): """ - cv.check_type('filters', filters, Iterable, openmc.FilterMeta) + cv.check_type('filters', filters, Iterable, FilterMeta) cv.check_type('filter_bins', filter_bins, Iterable, tuple) # If user did not specify any specific Filters, use them all @@ -1787,7 +1804,7 @@ def get_score_indices(self, scores): """ for score in scores: - if not isinstance(score, (str, openmc.CrossScore)): + if not isinstance(score, (str, CrossScore)): msg = f'Unable to get score indices for score "{score}" in ' \ f'ID="{self.id}" since it is not a string or CrossScore ' \ 'Tally' @@ -1984,9 +2001,9 @@ def get_pandas_dataframe(self, filters=True, nuclides=True, scores=True, column_name = 'score' for score in self.scores: - if isinstance(score, (str, openmc.CrossScore)): + if isinstance(score, (str, CrossScore)): scores.append(str(score)) - elif isinstance(score, openmc.AggregateScore): + elif isinstance(score, AggregateScore): scores.append(score.name) column_name = f'{score.aggregate_op}(score)' @@ -2086,7 +2103,7 @@ def get_reshaped_data(self, value='mean', expand_dims=False): for i, f in enumerate(self.filters): if expand_dims: # Mesh filter indices are backwards so we need to flip them - if type(f) in {openmc.MeshFilter, openmc.MeshBornFilter}: + if type(f) in {MeshFilter, MeshBornFilter}: fshape = f.shape[::-1] new_shape += fshape idx0, idx1 = i, i + len(fshape) - 1 @@ -2273,7 +2290,7 @@ def hybrid_product(self, other, binary_op, filter_product=None, else: all_filters = [self_copy.filters, other_copy.filters] for self_filter, other_filter in product(*all_filters): - new_filter = openmc.CrossFilter(self_filter, other_filter, + new_filter = CrossFilter(self_filter, other_filter, binary_op) new_tally.filters.append(new_filter) @@ -2284,7 +2301,7 @@ def hybrid_product(self, other, binary_op, filter_product=None, else: all_nuclides = [self_copy.nuclides, other_copy.nuclides] for self_nuclide, other_nuclide in product(*all_nuclides): - new_nuclide = openmc.CrossNuclide(self_nuclide, other_nuclide, + new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op) new_tally.nuclides.append(new_nuclide) @@ -2295,9 +2312,9 @@ def cross_score(score1, score2, binary_op): if score1 == score2: return score1 else: - return openmc.CrossScore(score1, score2, binary_op) + return CrossScore(score1, score2, binary_op) else: - return openmc.CrossScore(score1, score2, binary_op) + return CrossScore(score1, score2, binary_op) # Add scores to the new tally if score_product == 'entrywise': @@ -2506,16 +2523,16 @@ def _swap_filters(self, filter1, filter2): # Construct lists of tuples for the bins in each of the two filters filters = [type(filter1), type(filter2)] - if isinstance(filter1, openmc.DistribcellFilter): + if isinstance(filter1, DistribcellFilter): filter1_bins = [b for b in range(filter1.num_bins)] - elif isinstance(filter1, openmc.EnergyFunctionFilter): + elif isinstance(filter1, EnergyFunctionFilter): filter1_bins = [None] else: filter1_bins = filter1.bins - if isinstance(filter2, openmc.DistribcellFilter): + if isinstance(filter2, DistribcellFilter): filter2_bins = [b for b in range(filter2.num_bins)] - elif isinstance(filter2, openmc.EnergyFunctionFilter): + elif isinstance(filter2, EnergyFunctionFilter): filter2_bins = [None] else: filter2_bins = filter2.bins @@ -2648,11 +2665,11 @@ def _swap_scores(self, score1, score2): raise ValueError(msg) # Check that the scores are valid - if not isinstance(score1, (str, openmc.CrossScore)): + if not isinstance(score1, (str, CrossScore)): msg = 'Unable to swap score1 "{}" in Tally ID="{}" since it is ' \ 'not a string or CrossScore'.format(score1, self.id) raise ValueError(msg) - elif not isinstance(score2, (str, openmc.CrossScore)): + elif not isinstance(score2, (str, CrossScore)): msg = 'Unable to swap score2 "{}" in Tally ID="{}" since it is ' \ 'not a string or CrossScore'.format(score2, self.id) raise ValueError(msg) @@ -3296,7 +3313,7 @@ def get_slice(self, scores=[], filters=[], filter_bins=[], nuclides=[], new_filter.bins = [f.bins[i] for i in bin_indices] # Set number of bins manually for mesh/distribcell filters - if filter_type is openmc.DistribcellFilter: + if filter_type is DistribcellFilter: new_filter._num_bins = f._num_bins # Replace existing filter with new one @@ -3362,16 +3379,16 @@ def summation(self, scores=[], filter_type=None, std_dev = self.get_reshaped_data(value='std_dev') # Sum across any filter bins specified by the user - if isinstance(filter_type, openmc.FilterMeta): + if isinstance(filter_type, FilterMeta): find_filter = self.find_filter(filter_type) # If user did not specify filter bins, sum across all bins if len(filter_bins) == 0: bin_indices = np.arange(find_filter.num_bins) - if isinstance(find_filter, openmc.DistribcellFilter): + if isinstance(find_filter, DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) - elif isinstance(find_filter, openmc.EnergyFunctionFilter): + elif isinstance(find_filter, EnergyFunctionFilter): filter_bins = [None] else: filter_bins = find_filter.bins @@ -3400,7 +3417,7 @@ def summation(self, scores=[], filter_type=None, # Add AggregateFilter to the tally sum if not remove_filter: - filter_sum = openmc.AggregateFilter(self_filter, + filter_sum = AggregateFilter(self_filter, [tuple(filter_bins)], 'sum') tally_sum.filters.append(filter_sum) @@ -3423,7 +3440,7 @@ def summation(self, scores=[], filter_type=None, std_dev = np.sqrt(std_dev) # Add AggregateNuclide to the tally sum - nuclide_sum = openmc.AggregateNuclide(nuclides, 'sum') + nuclide_sum = AggregateNuclide(nuclides, 'sum') tally_sum.nuclides.append(nuclide_sum) # Add a copy of this tally's nuclides to the tally sum @@ -3441,7 +3458,7 @@ def summation(self, scores=[], filter_type=None, std_dev = np.sqrt(std_dev) # Add AggregateScore to the tally sum - score_sum = openmc.AggregateScore(scores, 'sum') + score_sum = AggregateScore(scores, 'sum') tally_sum.scores.append(score_sum) # Add a copy of this tally's scores to the tally sum @@ -3514,16 +3531,16 @@ def average(self, scores=[], filter_type=None, std_dev = self.get_reshaped_data(value='std_dev') # Average across any filter bins specified by the user - if isinstance(filter_type, openmc.FilterMeta): + if isinstance(filter_type, FilterMeta): find_filter = self.find_filter(filter_type) # If user did not specify filter bins, average across all bins if len(filter_bins) == 0: bin_indices = np.arange(find_filter.num_bins) - if isinstance(find_filter, openmc.DistribcellFilter): + if isinstance(find_filter, DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) - elif isinstance(find_filter, openmc.EnergyFunctionFilter): + elif isinstance(find_filter, EnergyFunctionFilter): filter_bins = [None] else: filter_bins = find_filter.bins @@ -3553,7 +3570,7 @@ def average(self, scores=[], filter_type=None, # Add AggregateFilter to the tally avg if not remove_filter: - filter_sum = openmc.AggregateFilter(self_filter, + filter_sum = AggregateFilter(self_filter, [tuple(filter_bins)], 'avg') tally_avg.filters.append(filter_sum) @@ -3577,7 +3594,7 @@ def average(self, scores=[], filter_type=None, std_dev = np.sqrt(std_dev) # Add AggregateNuclide to the tally avg - nuclide_avg = openmc.AggregateNuclide(nuclides, 'avg') + nuclide_avg = AggregateNuclide(nuclides, 'avg') tally_avg.nuclides.append(nuclide_avg) # Add a copy of this tally's nuclides to the tally avg @@ -3596,7 +3613,7 @@ def average(self, scores=[], filter_type=None, std_dev = np.sqrt(std_dev) # Add AggregateScore to the tally avg - score_sum = openmc.AggregateScore(scores, 'avg') + score_sum = AggregateScore(scores, 'avg') tally_avg.scores.append(score_sum) # Add a copy of this tally's scores to the tally avg @@ -3786,7 +3803,7 @@ def _create_mesh_subelements(self, root_element, memo=None): already_written = memo if memo else set() for tally in self: for f in tally.filters: - if isinstance(f, openmc.MeshFilter): + if isinstance(f, MeshFilter): if f.mesh.id in already_written: continue if len(f.mesh.name) > 0: @@ -3881,7 +3898,7 @@ def from_xml_element(cls, elem, meshes=None): # Read filter elements filters = {} for e in elem.findall('filter'): - filter = openmc.Filter.from_xml_element(e, meshes=meshes) + filter = Filter.from_xml_element(e, meshes=meshes) filters[filter.id] = filter # Read derivative elements diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 7797986df03..63af2596efc 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -11,6 +11,7 @@ import openmc from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh +from openmc.tallies import Tallies import openmc.checkvalue as cv from openmc.checkvalue import PathLike from ._xml import get_elem_list, get_text, clean_indentation @@ -499,6 +500,8 @@ class WeightWindowGenerator: Particle type the weight windows apply to method : {'magic', 'fw_cadis'} The weight window generation methodology applied during an update. + targets : :class:`openmc.Tallies` or iterable of int + Target tallies for local variance reduction via FW-CADIS. max_realizations : int The upper limit for number of tally realizations when generating weight windows. @@ -518,6 +521,8 @@ class WeightWindowGenerator: Particle type the weight windows apply to method : {'magic', 'fw_cadis'} The weight window generation methodology applied during an update. + targets : :class:`openmc.Tallies` or numpy.ndarray + Target tallies for local variance reduction via FW-CADIS. max_realizations : int The upper limit for number of tally realizations when generating weight windows. @@ -529,7 +534,7 @@ class WeightWindowGenerator: Whether or not to apply weight windows on the fly. """ - _MAGIC_PARAMS = {'value': str, 'threshold': float, 'ratio': float} + _WWG_PARAMS = {'value': str, 'threshold': float, 'ratio': float} def __init__( self, @@ -537,6 +542,7 @@ def __init__( energy_bounds: Sequence[float] | None = None, particle_type: str | int | openmc.ParticleType = 'neutron', method: str = 'magic', + targets: openmc.Tallies | Iterable[int] | None = None, max_realizations: int = 1, update_interval: int = 1, on_the_fly: bool = True @@ -549,6 +555,7 @@ def __init__( self.energy_bounds = energy_bounds self.particle_type = particle_type self.method = method + self.targets = targets self.max_realizations = max_realizations self.update_interval = update_interval self.on_the_fly = on_the_fly @@ -611,6 +618,22 @@ def method(self, m: str): self._check_update_parameters() except (TypeError, KeyError): warnings.warn(f'Update parameters are invalid for the "{m}" method.') + + @property + def targets(self) -> openmc.Tallies: + return self._targets + + @targets.setter + def targets(self, t): + if t is None: + self._targets = t + else: + cv.check_type('Local FW-CADIS target tallies', t, Iterable) + cv.check_greater_than('Local FW-CADIS target tallies', len(t), 0) + if not isinstance(t, openmc.Tallies): + cv.check_iterable_type('Local FW-CADIS target tallies', t, int) + t = np.asarray(list(t), dtype=int) + self._targets = t @property def max_realizations(self) -> int: @@ -638,13 +661,13 @@ def update_parameters(self) -> dict: def _check_update_parameters(self, params: dict): if self.method == 'magic' or self.method == 'fw_cadis': - check_params = self._MAGIC_PARAMS + check_params = self._WWG_PARAMS for key, val in params.items(): if key not in check_params: raise ValueError(f'Invalid param "{key}" for {self.method} ' 'weight window generation') - cv.check_type(f'weight window generation param: "{key}"', val, self._MAGIC_PARAMS[key]) + cv.check_type(f'weight window generation param: "{key}"', val, self._WWG_PARAMS[key]) @update_parameters.setter def update_parameters(self, params: dict): @@ -681,7 +704,7 @@ def _sanitize_update_parameters(cls, method: str, update_parameters: dict): The update parameters as-read from the XML node (keys: str, values: str) """ if method == 'magic' or method == 'fw_cadis': - check_params = cls._MAGIC_PARAMS + check_params = cls._WWG_PARAMS for param, param_type in check_params.items(): if param in update_parameters: @@ -707,6 +730,20 @@ def to_xml_element(self): otf_elem.text = str(self.on_the_fly).lower() method_elem = ET.SubElement(element, 'method') method_elem.text = self.method + if self.targets is not None: + if self.method != 'fw_cadis': + raise ValueError( + "FW-CADIS update method is required in order to use " \ + "target tallies for WeightWindowGenerator.") + elif isinstance(self.targets, openmc.Tallies): + raise RuntimeError( + "FW-CADIS target tallies must be checked to ensure they are " \ + "present on model.tallies. Use model.export_to_xml() or " \ + "model.export_to_model_xml() to link FW-CADIS target tallies.") + else: + targets_elem = ET.SubElement(element, 'targets') + targets_elem.text = ' '.join(str(tally_id) for tally_id in self.targets) + if self.update_parameters is not None: self._update_parameters_subelement(element) @@ -733,8 +770,8 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self: mesh_id = int(get_text(elem, 'mesh')) mesh = meshes[mesh_id] - - energy_bounds = get_elem_list(elem, "energy_bounds, float") + + energy_bounds = get_elem_list(elem, "energy_bounds", float) particle_type = get_text(elem, 'particle_type') wwg = cls(mesh, energy_bounds, particle_type) @@ -743,6 +780,14 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self: wwg.update_interval = int(get_text(elem, 'update_interval')) wwg.on_the_fly = bool(get_text(elem, 'on_the_fly')) wwg.method = get_text(elem, 'method') + targets_elem = elem.find('targets') + if targets_elem is not None: + if wwg.method != 'fw_cadis': + raise ValueError( + "FW-CADIS update method is required in order to use " \ + "target tallies for WeightWindowGenerator.") + else: + wwg.targets = get_elem_list(elem, "targets") if elem.find('update_parameters') is not None: update_parameters = {} diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 1a6e7c0be46..06c6ef14d71 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -31,9 +31,11 @@ RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ { RandomRayVolumeEstimator::HYBRID}; bool FlatSourceDomain::volume_normalized_flux_tallies_ {false}; bool FlatSourceDomain::adjoint_ {false}; +bool FlatSourceDomain::fw_cadis_local_ {false}; double FlatSourceDomain::diagonal_stabilization_rho_ {1.0}; std::unordered_map>> FlatSourceDomain::mesh_domain_map_; +std::vector FlatSourceDomain::fw_cadis_local_targets_; FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) { @@ -1000,7 +1002,9 @@ void FlatSourceDomain::output_to_vtk() const void FlatSourceDomain::apply_external_source_to_source_region( int src_idx, SourceRegionHandle& srh) { - auto s = model::external_sources[src_idx].get(); + auto s = (adjoint_ && !model::adjoint_sources.empty()) + ? model::adjoint_sources[src_idx].get() + : model::external_sources[src_idx].get(); auto is = dynamic_cast(s); auto discrete = dynamic_cast(is->energy()); double strength_factor = is->strength(); @@ -1071,13 +1075,17 @@ void FlatSourceDomain::count_external_source_regions() } } -void FlatSourceDomain::convert_external_sources() +void FlatSourceDomain::convert_external_sources(bool use_adjoint_sources) { + // Determine whether forward or (local) adjoint sources are desired + const auto& sources = + use_adjoint_sources ? model::adjoint_sources : model::external_sources; + // Loop over external sources - for (int es = 0; es < model::external_sources.size(); es++) { + for (int es = 0; es < sources.size(); es++) { // Extract source information - Source* s = model::external_sources[es].get(); + Source* s = sources[es].get(); IndependentSource* is = dynamic_cast(s); Discrete* energy = dynamic_cast(is->energy()); const std::unordered_set& domain_ids = is->domain_ids(); @@ -1223,7 +1231,7 @@ void FlatSourceDomain::flatten_xs() } } -void FlatSourceDomain::set_adjoint_sources() +void FlatSourceDomain::set_fw_adjoint_sources() { // Set the adjoint external source to 1/forward_flux. If the forward flux is // negative, zero, or extremely close to zero, set the adjoint source to zero, @@ -1252,6 +1260,10 @@ void FlatSourceDomain::set_adjoint_sources() source_regions_.external_source(sr, g) = 0.0; } else { source_regions_.external_source(sr, g) = 1.0 / flux; + if (!std::isfinite(source_regions_.external_source(sr, g))) { + // If the flux is NaN or Inf, set the adjoint source to zero + source_regions_.external_source(sr, g) = 0.0; + } } if (flux > 0.0) { source_regions_.external_source_present(sr) = 1; @@ -1283,6 +1295,7 @@ void FlatSourceDomain::set_adjoint_sources() source_regions_.external_source_present(sr) = 0; } } + // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for @@ -1297,8 +1310,86 @@ void FlatSourceDomain::set_adjoint_sources() sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] * source_regions_.density_mult(sr); source_regions_.external_source(sr, g) /= sigma_t; + if (!std::isfinite(source_regions_.external_source(sr, g))) { + // If the flux is NaN or Inf, set the adjoint source to zero + source_regions_.external_source(sr, g) = 0.0; + } } } + + if (fw_cadis_local_) { +// Only external sources that have a non-mesh type tally task should remain +// non-zero. Everything else gets zero'd out. +#pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions(); sr++) { + + // If there is already no external source, don't need to do anything + if (source_regions_.external_source_present(sr) == 0) { + continue; + } + + // If there is an adjoint source term here, then we need to check it. + + // We will track if ANY group has a valid local FW-CADIS source term + bool has_any_sources = false; + + // Now, loop over groups + for (int g = 0; g < negroups_; g++) { + + // If there are no tally tasks associated with this source element + // then it is not a local FW-CADIS source, so we continue to the next + // group + if (source_regions_.tally_task(sr, g).empty()) { + source_regions_.external_source(sr, g) = 0.0; + continue; + } + + // If there are tally tasks, we can through them and check if + // any of them are local FW-CADIS targets. + + // We track if ANY of the tasks are local FW-CADIS target tallies + bool local_fw_cadis_target_region = false; + + // Now we loop through + for (const auto& task : source_regions_.tally_task(sr, g)) { + Tally& tally {*model::tallies[task.tally_idx]}; + const auto t_id = tally.id(); + + // Search for target tallies + if (std::find(fw_cadis_local_targets_.begin(), + fw_cadis_local_targets_.end(), + t_id) != fw_cadis_local_targets_.end()) { + local_fw_cadis_target_region = true; + break; + } + } + + // If ANY of the tasks is a local FW-CADIS target, + // Then we keep the source term and set that this + // source region has a valid FW-CADIS source term. + // Otherwise, we zero out the source term. + if (local_fw_cadis_target_region) { + has_any_sources = true; + } else { + source_regions_.external_source(sr, g) = 0.0; + } + } // End loop over groups + + // If there were any valid FW-CADIS source terms for any + // of the groups, then the SR as a whole counts as a source + if (has_any_sources) { + source_regions_.external_source_present(sr) = 1; + } else { + source_regions_.external_source_present(sr) = 0; + } + } // End loop over source regions + } // End local FW-CADIS logic +} + +void FlatSourceDomain::set_local_adjoint_sources() +{ + // Set the external source to user-specified adjoint sources. + convert_external_sources(true); } void FlatSourceDomain::transpose_scattering_matrix() diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 80cbfc3fe5b..24650afb8dd 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -168,14 +168,12 @@ void validate_random_ray_inputs() "constrained by domain id (cell, material, or universe) in " "random ray mode."); } else if (is->domain_ids().size() > 0 && sp) { - // If both a domain constraint and a non-default point source location - // are specified, notify user that domain constraint takes precedence. - if (sp->r().x == 0.0 && sp->r().y == 0.0 && sp->r().z == 0.0) { - warning("Fixed source has both a domain constraint and a point " - "type spatial distribution. The domain constraint takes " - "precedence in random ray mode -- point source coordinate " - "will be ignored."); - } + // If both a domain constraint and a point source location are + // specified, notify user that domain constraint takes precedence. + warning("Fixed source has both a domain constraint and a point " + "type spatial distribution. The domain constraint takes " + "precedence in random ray mode -- point source coordinate " + "will be ignored."); } // Check that a discrete energy distribution was used @@ -189,6 +187,56 @@ void validate_random_ray_inputs() } } + // Validate adjoint sources + /////////////////////////////////////////////////////////////////// + if (FlatSourceDomain::adjoint_ && !model::adjoint_sources.empty()) { + for (int i = 0; i < model::adjoint_sources.size(); i++) { + Source* s = model::adjoint_sources[i].get(); + + // Check for independent source + IndependentSource* is = dynamic_cast(s); + + if (!is) { + fatal_error( + "Only IndependentSource adjoint source types are allowed in " + "random ray mode"); + } + + // Check for isotropic source + UnitSphereDistribution* angle_dist = is->angle(); + Isotropic* id = dynamic_cast(angle_dist); + if (!id) { + fatal_error( + "Invalid source definition -- only isotropic adjoint sources are " + "allowed in random ray mode."); + } + + // Validate that a domain ID was specified OR that it is a point source + auto sp = dynamic_cast(is->space()); + if (is->domain_ids().size() == 0 && !sp) { + fatal_error("Adjoint sources must be point source or spatially " + "constrained by domain id (cell, material, or universe) in " + "random ray mode."); + } else if (is->domain_ids().size() > 0 && sp) { + // If both a domain constraint and a point source location are + // specified, notify user that domain constraint takes precedence. + warning("Adjoint source has both a domain constraint and a point " + "type spatial distribution. The domain constraint takes " + "precedence in random ray mode -- point source coordinate " + "will be ignored."); + } + + // Check that a discrete energy distribution was used + Distribution* d = is->energy(); + Discrete* dd = dynamic_cast(d); + if (!dd) { + fatal_error( + "Only discrete (multigroup) energy distributions are allowed for " + "adjoint sources in random ray mode."); + } + } + } + // Validate plotting files /////////////////////////////////////////////////////////////////// for (int p = 0; p < model::plots.size(); p++) { @@ -232,20 +280,22 @@ void validate_random_ray_inputs() warning( "Linear sources may result in negative fluxes in small source regions " "generated by mesh subdivision. Negative sources may result in low " - "quality FW-CADIS weight windows. We recommend you use flat source mode " - "when generating weight windows with an overlaid mesh tally."); + "quality FW-CADIS weight windows. We recommend you use flat source " + "mode when generating weight windows with an overlaid mesh tally."); } } -void print_adjoint_header() +void openmc_finalize_random_ray() { - if (!FlatSourceDomain::adjoint_) - // If we're going to do an adjoint simulation afterwards, report that this - // is the initial forward flux solve. - header("FORWARD FLUX SOLVE", 3); - else - // Otherwise report that we are doing the adjoint simulation - header("ADJOINT FLUX SOLVE", 3); + FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID; + FlatSourceDomain::volume_normalized_flux_tallies_ = false; + FlatSourceDomain::adjoint_ = false; + FlatSourceDomain::fw_cadis_local_ = false; + FlatSourceDomain::fw_cadis_local_targets_.clear(); + FlatSourceDomain::mesh_domain_map_.clear(); + RandomRay::ray_source_.reset(); + RandomRay::source_shape_ = RandomRaySourceShape::FLAT; + RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } //============================================================================== @@ -278,16 +328,6 @@ RandomRaySimulation::RandomRaySimulation() // Convert OpenMC native MGXS into a more efficient format // internal to the random ray solver domain_->flatten_xs(); - - // Check if adjoint calculation is needed. If it is, we will run the forward - // calculation first and then the adjoint calculation later. - adjoint_needed_ = FlatSourceDomain::adjoint_; - - // Adjoint is always false for the forward calculation - FlatSourceDomain::adjoint_ = false; - - // The first simulation is run after initialization - is_first_simulation_ = true; } void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() @@ -295,29 +335,51 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() domain_->apply_meshes(); if (settings::run_mode == RunMode::FIXED_SOURCE) { // Transfer external source user inputs onto random ray source regions - domain_->convert_external_sources(); + domain_->convert_external_sources(false); domain_->count_external_source_regions(); } } -void RandomRaySimulation::prepare_fixed_sources_adjoint() +void RandomRaySimulation::prepare_fw_fixed_sources_adjoint() { + // Prepare adjoint fixed sources using forward flux domain_->source_regions_.adjoint_reset(); if (settings::run_mode == RunMode::FIXED_SOURCE) { - domain_->set_adjoint_sources(); + domain_->set_fw_adjoint_sources(); + } +} + +void RandomRaySimulation::prepare_local_fixed_sources_adjoint() +{ + if (settings::run_mode == RunMode::FIXED_SOURCE) { + domain_->set_local_adjoint_sources(); } } -void RandomRaySimulation::prepare_adjoint_simulation() +void RandomRaySimulation::prepare_adjoint_simulation(bool fw_adjoint) { - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; + reset_timers(); - // Reset k-eff - domain_->k_eff_ = 1.0; + if (mpi::master) + header("ADJOINT FLUX SOLVE", 3); - // Initialize adjoint fixed sources, if present - prepare_fixed_sources_adjoint(); + if (fw_adjoint) { + // Forward simulation has already been run; + // Configure the domain for adjoint simulation and + // re-initialize OpenMC general data structures + FlatSourceDomain::adjoint_ = true; + + openmc_simulation_init(); + + prepare_fw_fixed_sources_adjoint(); + } else { + // Initialize adjoint fixed sources + domain_->apply_meshes(); + prepare_local_fixed_sources_adjoint(); + domain_->count_external_source_regions(); + } + + domain_->k_eff_ = 1.0; // Transpose scattering matrix domain_->transpose_scattering_matrix(); @@ -328,18 +390,6 @@ void RandomRaySimulation::prepare_adjoint_simulation() void RandomRaySimulation::simulate() { - if (!is_first_simulation_) { - if (mpi::master && adjoint_needed_) - openmc::print_adjoint_header(); - - // Reset the timers and reinitialize the general OpenMC datastructures if - // this is after the first simulation - reset_timers(); - - // Initialize OpenMC general data structures - openmc_simulation_init(); - } - // Begin main simulation timer simulation::time_total.start(); @@ -435,7 +485,7 @@ void RandomRaySimulation::simulate() // End main simulation timer simulation::time_total.stop(); - // Normalize and save the final flux + // Normalize and save the final forward flux double source_normalization_factor = domain_->compute_fixed_source_normalization_factor() / (settings::n_batches - settings::n_inactive); @@ -451,11 +501,6 @@ void RandomRaySimulation::simulate() // Output all simulation results output_simulation_results(); - - // Toggle that the simulation object has been initialized after the first - // simulation - if (is_first_simulation_) - is_first_simulation_ = false; } void RandomRaySimulation::output_simulation_results() const @@ -622,17 +667,6 @@ void RandomRaySimulation::print_results_random_ray( } } -void openmc_finalize_random_ray() -{ - FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID; - FlatSourceDomain::volume_normalized_flux_tallies_ = false; - FlatSourceDomain::adjoint_ = false; - FlatSourceDomain::mesh_domain_map_.clear(); - RandomRay::ray_source_.reset(); - RandomRay::source_shape_ = RandomRaySourceShape::FLAT; - RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; -} - } // namespace openmc //============================================================================== @@ -645,12 +679,25 @@ void openmc_run_random_ray() // Run forward simulation ////////////////////////////////////////////////////////// - if (openmc::mpi::master) { - if (openmc::FlatSourceDomain::adjoint_) { - openmc::FlatSourceDomain::adjoint_ = false; - openmc::print_adjoint_header(); - openmc::FlatSourceDomain::adjoint_ = true; - } + // Check if adjoint calculation is needed, and if local adjoint source(s) + // are present. If an adjoint calculation is needed and no sources are + // specified, we will run a forward calculation first to calculate adjoint + // sources for global variance reduction, then perform an adjoint + // calculation later. + bool adjoint_needed = openmc::FlatSourceDomain::adjoint_; + bool fw_adjoint = openmc::model::adjoint_sources.empty() && adjoint_needed; + + // If we're going to do an adjoint simulation with forward-weighted adjoint + // sources afterwards, report that this is the initial forward flux solve. + if (!adjoint_needed || fw_adjoint) { + // Configure the domain for forward simulation + openmc::FlatSourceDomain::adjoint_ = false; + + if (adjoint_needed && openmc::mpi::master) + openmc::header("FORWARD FLUX SOLVE", 3); + } else { + // Configure domain for adjoint simulation (later) + openmc::FlatSourceDomain::adjoint_ = true; } // Initialize OpenMC general data structures @@ -663,21 +710,25 @@ void openmc_run_random_ray() // Initialize Random Ray Simulation Object openmc::RandomRaySimulation sim; - // Initialize fixed sources, if present - sim.apply_fixed_sources_and_mesh_domains(); + if (!adjoint_needed || fw_adjoint) { + // Initialize fixed sources, if present + sim.apply_fixed_sources_and_mesh_domains(); - // Run initial random ray simulation - sim.simulate(); + // Execute random ray simulation + sim.simulate(); + } ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (sim.adjoint_needed_) { - // Setup for adjoint simulation - sim.prepare_adjoint_simulation(); - - // Run adjoint simulation - sim.simulate(); + if (!adjoint_needed) { + return; } + + // Setup for adjoint simulation + sim.prepare_adjoint_simulation(fw_adjoint); + + // Execute random ray simulation + sim.simulate(); } diff --git a/src/settings.cpp b/src/settings.cpp index ab9f9a5aafa..bb6991da955 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -280,8 +280,9 @@ void get_run_parameters(pugi::xml_node node_base) } else { fatal_error("Specify random ray inactive distance in settings XML"); } - if (check_for_node(random_ray_node, "source")) { - xml_node source_node = random_ray_node.child("source"); + if (check_for_node(random_ray_node, "ray_source")) { + xml_node ray_source_node = random_ray_node.child("ray_source"); + xml_node source_node = ray_source_node.child("source"); // Get point to list of elements and make sure there is at least // one RandomRay::ray_source_ = Source::create(source_node); @@ -369,6 +370,13 @@ void get_run_parameters(pugi::xml_node node_base) "between 0 and 1"); } } + if (check_for_node(random_ray_node, "adjoint_source")) { + pugi::xml_node adj_source_node = random_ray_node.child("adjoint_source"); + for (pugi::xml_node source_node : adj_source_node.children("source")) { + // Find any local adjoint sources + model::adjoint_sources.push_back(Source::create(source_node)); + } + } } } @@ -1276,6 +1284,16 @@ void read_settings_xml(pugi::xml_node root) break; } } + // If any weight window generators have local FW-CADIS target tallies, + // user-defined adjoint sources cannot be used at the same time. + if (!model::adjoint_sources.empty()) { + for (const auto& wwg : variance_reduction::weight_windows_generators) { + if (!wwg->targets_.empty()) { + fatal_error("Cannot use both user-defined adjoint sources and " + "FW-CADIS target tallies at the same time."); + } + } + } } // Set up weight window checkpoints diff --git a/src/source.cpp b/src/source.cpp index f0b8f48edda..524a72bf0e4 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -62,6 +62,8 @@ namespace model { vector> external_sources; +vector> adjoint_sources; + DiscreteIndex external_sources_probability; } // namespace model @@ -710,6 +712,7 @@ SourceSite sample_external_source(uint64_t* seed) void free_memory_source() { model::external_sources.clear(); + model::adjoint_sources.clear(); reset_source_rejection_counters(); } diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index c00872e5643..0614110cd32 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -621,7 +621,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, } } } else { - // For FW-CADIS, weight windows are inversely proportional to the adjoint + // For (FW-)CADIS, weight windows are inversely proportional to the adjoint // fluxes. We normalize the weight windows across all energy groups. #pragma omp parallel for collapse(2) schedule(static) for (int e = 0; e < e_bins; e++) { @@ -801,6 +801,13 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) fatal_error("FW-CADIS can only be run in random ray solver mode."); } FlatSourceDomain::adjoint_ = true; + if (check_for_node(node, "targets")) { + FlatSourceDomain::fw_cadis_local_ = true; + targets_ = get_node_array(node, "targets"); + FlatSourceDomain::fw_cadis_local_targets_.insert( + std::end(FlatSourceDomain::fw_cadis_local_targets_), + std::begin(targets_), std::end(targets_)); + } } else { fatal_error(fmt::format( "Unknown weight window update method '{}' specified", method_string)); diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat index 0adfc548848..94e2709768f 100644 --- a/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat +++ b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true true naive diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat index 073348c41e7..755afd6c439 100644 --- a/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat +++ b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true true diff --git a/tests/regression_tests/random_ray_adjoint_local/__init__.py b/tests/regression_tests/random_ray_adjoint_local/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_adjoint_local/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_local/inputs_true.dat new file mode 100644 index 00000000000..9021d1675d1 --- /dev/null +++ b/tests/regression_tests/random_ray_adjoint_local/inputs_true.dat @@ -0,0 +1,293 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + + + + + + fixed source + 500 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 800.0 + 100.0 + + + + 0.0 0.0 0.0 35.0 35.0 35.0 + + + + true + + + + + + true + + + + 100.0 1.0 + + + cell + 6 7 + + + + naive + + + 14 14 14 + 0.0 0.0 0.0 + 35.0 35.0 35.0 + + + + + 6 + + + 7 + + + 3 + + + 2 + + + 1 + + + 1 + flux + tracklength + + + 2 + flux + tracklength + + + 3 + flux + tracklength + + + 4 + flux + tracklength + + + 5 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_adjoint_local/results_true.dat b/tests/regression_tests/random_ray_adjoint_local/results_true.dat new file mode 100644 index 00000000000..daa94856559 --- /dev/null +++ b/tests/regression_tests/random_ray_adjoint_local/results_true.dat @@ -0,0 +1,15 @@ +tally 1: +2.215273E+01 +9.815738E+01 +tally 2: +1.873933E+01 +7.023420E+01 +tally 3: +4.802282E-01 +4.612707E-02 +tally 4: +2.516720E-01 +1.271063E-02 +tally 5: +1.169938E-02 +3.277334E-05 diff --git a/tests/regression_tests/random_ray_adjoint_local/test.py b/tests/regression_tests/random_ray_adjoint_local/test.py new file mode 100644 index 00000000000..c11b8e84700 --- /dev/null +++ b/tests/regression_tests/random_ray_adjoint_local/test.py @@ -0,0 +1,35 @@ +import os +import openmc + +from openmc.examples import random_ray_three_region_cube_with_detectors + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_adjoint_local(): + model = random_ray_three_region_cube_with_detectors() + + detector1_cells = model.geometry.get_cells_by_name("detector 1") + detector2_cells = model.geometry.get_cells_by_name("detector 2") + detector_cells = detector1_cells + detector2_cells + + strengths = [1.0] + midpoints = [100.0] + energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths) + + adj_source = openmc.IndependentSource(energy=energy_distribution, constraints={ + 'domains': detector_cells}, strength=3.14) + + model.settings.random_ray['adjoint'] = True + model.settings.random_ray['adjoint_source'] = adj_source + model.settings.random_ray['volume_estimator'] = 'naive' + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat index 464c89a5df9..86d5ec4abd5 100644 --- a/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat index 464c89a5df9..86d5ec4abd5 100644 --- a/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat index 464c89a5df9..86d5ec4abd5 100644 --- a/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat index 9f3a827f687..b00935ef38a 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat index edf68f7e25b..472406fa882 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat index edf68f7e25b..472406fa882 100644 --- a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat index 80a166c6786..15981f7fa5f 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/model/inputs_true.dat @@ -38,11 +38,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat index 464c89a5df9..86d5ec4abd5 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/infinite_medium/user/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat index 80a166c6786..15981f7fa5f 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/model/inputs_true.dat @@ -38,11 +38,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat index 464c89a5df9..86d5ec4abd5 100644 --- a/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_source_energy/stochastic_slab/user/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat index 08a30176bef..c60e6a04199 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/infinite_medium/inputs_true.dat @@ -47,11 +47,13 @@ 200.0 400.0 200.0 - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat index 08a30176bef..c60e6a04199 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/material_wise/inputs_true.dat @@ -47,11 +47,13 @@ 200.0 400.0 200.0 - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat index 08a30176bef..c60e6a04199 100644 --- a/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_temperature/stochastic_slab/inputs_true.dat @@ -47,11 +47,13 @@ 200.0 400.0 200.0 - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat index bcf4d0d90cb..eacd54f8354 100644 --- a/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat +++ b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat @@ -86,11 +86,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true diff --git a/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat index e90f25973e2..f369bae89f3 100644 --- a/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat +++ b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_cell_temperature/inputs_true.dat b/tests/regression_tests/random_ray_cell_temperature/inputs_true.dat index e8674e6eafb..99363ed8760 100644 --- a/tests/regression_tests/random_ray_cell_temperature/inputs_true.dat +++ b/tests/regression_tests/random_ray_cell_temperature/inputs_true.dat @@ -89,11 +89,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true diff --git a/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat b/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat index 47325ebd7d5..11100e88e12 100644 --- a/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat +++ b/tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat @@ -41,11 +41,13 @@ multi-group - - - -0.63 -0.63 -1.0 0.63 0.63 1.0 - - + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 150.0 diff --git a/tests/regression_tests/random_ray_entropy/settings.xml b/tests/regression_tests/random_ray_entropy/settings.xml index 81deaa7751d..0d830417b62 100644 --- a/tests/regression_tests/random_ray_entropy/settings.xml +++ b/tests/regression_tests/random_ray_entropy/settings.xml @@ -6,11 +6,13 @@ 5 multi-group - - - 0.0 0.0 0.0 100.0 100.0 100.0 - - + + + + 0.0 0.0 0.0 100.0 100.0 100.0 + + + 40.0 400.0 diff --git a/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat index 9f1987f3acc..d650bbaf95c 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat index b4f57dbfa8a..98a51add1f3 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat index ab91f74e50d..20deba664bd 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat index 220fa7db643..2268d82c391 100644 --- a/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat index f8c4430852f..fe95baa7bb5 100644 --- a/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear_xy diff --git a/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat index c84e544fcc4..a5632ece960 100644 --- a/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat index 05c4846e6b4..9d22603c63c 100644 --- a/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat index 0c870e10067..de941f10fbb 100644 --- a/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + false diff --git a/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat index ab91f74e50d..20deba664bd 100644 --- a/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat index 0c05a71df32..943468a1095 100644 --- a/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat @@ -110,11 +110,13 @@ 40.0 40.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + false flat diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat index a67495bf169..650953c4b06 100644 --- a/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat @@ -110,11 +110,13 @@ 40.0 40.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + false linear_xy diff --git a/tests/regression_tests/random_ray_halton_samples/inputs_true.dat b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat index 36d5f6f2276..1b86d2daeea 100644 --- a/tests/regression_tests/random_ray_halton_samples/inputs_true.dat +++ b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true halton diff --git a/tests/regression_tests/random_ray_k_eff/inputs_true.dat b/tests/regression_tests/random_ray_k_eff/inputs_true.dat index 545bd1d4578..72b783344fd 100644 --- a/tests/regression_tests/random_ray_k_eff/inputs_true.dat +++ b/tests/regression_tests/random_ray_k_eff/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true diff --git a/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat b/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat index 98badea18d8..f6e9c8e3e71 100644 --- a/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat +++ b/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true diff --git a/tests/regression_tests/random_ray_linear/linear/inputs_true.dat b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat index a43a66e71c7..269d9892ebf 100644 --- a/tests/regression_tests/random_ray_linear/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true linear diff --git a/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat index 7f76f2fd1c1..217e9551614 100644 --- a/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat +++ b/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat @@ -80,11 +80,13 @@ 100.0 20.0 - - - -1.26 -1.26 -1 1.26 1.26 1 - - + + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true linear_xy diff --git a/tests/regression_tests/random_ray_low_density/inputs_true.dat b/tests/regression_tests/random_ray_low_density/inputs_true.dat index ab91f74e50d..20deba664bd 100644 --- a/tests/regression_tests/random_ray_low_density/inputs_true.dat +++ b/tests/regression_tests/random_ray_low_density/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat index 088f803bfa8..b4bd263f5ac 100644 --- a/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat +++ b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat @@ -206,11 +206,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/random_ray_s2/inputs_true.dat b/tests/regression_tests/random_ray_s2/inputs_true.dat index aad8ea28b8c..c0dc6292f7f 100644 --- a/tests/regression_tests/random_ray_s2/inputs_true.dat +++ b/tests/regression_tests/random_ray_s2/inputs_true.dat @@ -37,11 +37,13 @@ 100.0 400.0 - - - 0.0 -5.0 -5.0 40.0 5.0 5.0 - - + + + + 0.0 -5.0 -5.0 40.0 5.0 5.0 + + + flat s2 diff --git a/tests/regression_tests/random_ray_void/flat/inputs_true.dat b/tests/regression_tests/random_ray_void/flat/inputs_true.dat index aa28e7b68bf..66390c76661 100644 --- a/tests/regression_tests/random_ray_void/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_void/flat/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true flat diff --git a/tests/regression_tests/random_ray_void/linear/inputs_true.dat b/tests/regression_tests/random_ray_void/linear/inputs_true.dat index e4b2f22fa27..45228a03955 100644 --- a/tests/regression_tests/random_ray_void/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_void/linear/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear diff --git a/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat index 8e8a8ed9b81..4d1af46b121 100644 --- a/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true hybrid diff --git a/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat index 1e25b97da66..a268d55d04a 100644 --- a/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true naive diff --git a/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat index 78c16269763..777ccaea510 100644 --- a/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true simulation_averaged diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat index 47a8a718249..dd11567f69d 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear hybrid diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat index 80a9ada4d5b..6933fba435e 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear naive diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat index 4f032a62a82..3ccab1d21b7 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat @@ -207,11 +207,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true linear simulation_averaged diff --git a/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat index 5fa6505ddf4..6bdabfbee31 100644 --- a/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat @@ -222,11 +222,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true naive diff --git a/tests/regression_tests/weightwindows_fw_cadis_local/__init__.py b/tests/regression_tests/weightwindows_fw_cadis_local/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/weightwindows_fw_cadis_local/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis_local/inputs_true.dat new file mode 100644 index 00000000000..ffdd977ec06 --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis_local/inputs_true.dat @@ -0,0 +1,273 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + + + + + + fixed source + 500 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + + 2 + neutron + 10 + 1 + true + fw_cadis + 1 2 + + + + 7 7 7 + 0.0 0.0 0.0 + 35.0 35.0 35.0 + + + 800.0 + 100.0 + + + + 0.0 0.0 0.0 35.0 35.0 35.0 + + + + true + + + + + + naive + + + 14 14 14 + 0.0 0.0 0.0 + 35.0 35.0 35.0 + + + + + 6 + + + 7 + + + 1 + flux + tracklength + + + 2 + flux + tracklength + + + diff --git a/tests/regression_tests/weightwindows_fw_cadis_local/results_true.dat b/tests/regression_tests/weightwindows_fw_cadis_local/results_true.dat new file mode 100644 index 00000000000..5991fc621c3 --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis_local/results_true.dat @@ -0,0 +1,696 @@ +RegularMesh + ID = 2 + Name = + Dimensions = 3 + Voxels = [7 7 7] + Lower left = [0. 0. 0.] + Upper Right = [np.float64(35.0), np.float64(35.0), np.float64(35.0)] + Width = [5. 5. 5.] +Lower Bounds +1.52e-01 +1.53e-01 +1.78e-01 +1.97e-01 +2.41e-01 +3.94e-01 +2.26e-03 +1.43e-01 +1.58e-01 +1.75e-01 +1.86e-01 +2.06e-01 +4.83e-02 +2.11e-03 +1.31e-01 +1.37e-01 +1.61e-01 +1.76e-01 +2.04e-01 +1.15e-01 +5.04e-03 +1.05e-01 +1.31e-01 +1.49e-01 +1.72e-01 +1.81e-01 +4.24e-02 +5.31e-03 +6.79e-02 +9.98e-02 +1.56e-01 +1.68e-01 +1.78e-01 +1.39e-02 +5.49e-03 +2.69e-03 +6.66e-03 +2.36e-02 +1.65e-02 +1.71e-02 +1.12e-02 +1.94e-03 +1.12e-04 +5.43e-04 +1.23e-03 +1.51e-03 +1.51e-03 +1.31e-03 +1.21e-03 +1.52e-01 +1.60e-01 +1.81e-01 +2.03e-01 +2.16e-01 +2.58e-02 +2.07e-03 +1.45e-01 +1.59e-01 +1.84e-01 +2.07e-01 +2.14e-01 +8.82e-02 +2.92e-03 +1.36e-01 +1.53e-01 +1.61e-01 +1.92e-01 +2.28e-01 +1.82e-01 +2.70e-03 +1.18e-01 +1.42e-01 +1.52e-01 +1.81e-01 +1.98e-01 +2.87e-02 +6.46e-03 +1.00e-01 +1.41e-01 +1.80e-01 +1.73e-01 +2.05e-01 +1.09e-01 +3.29e-03 +5.37e-03 +1.44e-02 +4.53e-02 +2.16e-02 +7.53e-02 +2.38e-02 +1.38e-03 +4.52e-04 +1.11e-03 +1.39e-03 +1.33e-03 +1.26e-03 +1.22e-03 +9.78e-04 +1.83e-01 +1.77e-01 +1.83e-01 +2.03e-01 +2.27e-01 +7.70e-02 +3.55e-03 +1.63e-01 +1.80e-01 +1.87e-01 +1.86e-01 +2.20e-01 +4.43e-02 +2.00e-03 +1.68e-01 +1.60e-01 +1.63e-01 +1.78e-01 +2.12e-01 +9.86e-02 +2.96e-03 +1.56e-01 +1.59e-01 +1.63e-01 +1.72e-01 +1.81e-01 +7.15e-02 +1.77e-03 +1.58e-01 +1.80e-01 +1.87e-01 +1.66e-01 +1.66e-01 +1.22e-02 +2.65e-03 +2.11e-02 +3.01e-02 +4.85e-02 +1.48e-02 +4.67e-02 +5.38e-03 +1.01e-03 +8.10e-04 +8.55e-04 +1.43e-03 +1.63e-03 +1.05e-03 +1.04e-03 +6.84e-04 +2.01e-01 +2.17e-01 +2.26e-01 +2.37e-01 +2.45e-01 +3.25e-02 +2.68e-02 +2.09e-01 +2.09e-01 +2.07e-01 +2.02e-01 +2.15e-01 +2.42e-02 +2.97e-03 +2.01e-01 +1.94e-01 +1.93e-01 +1.80e-01 +1.72e-01 +1.32e-02 +1.54e-03 +1.84e-01 +1.87e-01 +1.63e-01 +1.49e-01 +1.50e-01 +6.42e-02 +1.91e-03 +1.93e-01 +2.15e-01 +1.93e-01 +1.42e-01 +1.22e-01 +7.21e-03 +9.37e-04 +1.67e-02 +5.33e-02 +5.24e-02 +1.13e-02 +1.23e-02 +4.44e-03 +6.14e-04 +7.25e-04 +6.85e-04 +8.30e-04 +9.75e-04 +6.78e-04 +5.52e-04 +5.29e-04 +2.82e-01 +2.73e-01 +2.85e-01 +2.57e-01 +2.70e-01 +2.24e-02 +2.57e-03 +2.52e-01 +2.49e-01 +2.38e-01 +2.53e-01 +2.37e-01 +1.48e-02 +2.30e-03 +2.18e-01 +2.22e-01 +2.20e-01 +1.96e-01 +1.80e-01 +1.93e-02 +1.53e-03 +2.21e-01 +2.39e-01 +1.99e-01 +1.67e-01 +1.45e-01 +1.11e-02 +8.58e-04 +2.57e-01 +2.32e-01 +1.89e-01 +1.26e-01 +8.31e-02 +8.12e-03 +6.25e-04 +8.53e-02 +7.87e-02 +1.53e-02 +1.22e-02 +7.53e-03 +1.75e-03 +3.39e-04 +1.08e-03 +8.97e-04 +6.11e-04 +5.24e-04 +4.33e-04 +2.49e-04 +2.36e-04 +4.78e-01 +1.38e-01 +5.00e-01 +2.99e-02 +3.12e-02 +1.07e-01 +1.42e-03 +3.01e-01 +4.37e-02 +1.79e-01 +2.39e-01 +7.54e-02 +1.13e-02 +1.17e-03 +3.44e-01 +2.30e-01 +5.55e-02 +5.69e-02 +1.70e-02 +8.02e-03 +8.18e-04 +1.77e-01 +6.57e-02 +6.14e-02 +3.20e-02 +1.32e-02 +6.75e-03 +6.17e-04 +3.49e-02 +1.65e-02 +1.95e-02 +8.71e-03 +9.02e-03 +2.07e-03 +3.41e-04 +7.97e-02 +2.60e-02 +1.65e-02 +6.67e-03 +2.61e-03 +5.74e-04 +1.41e-04 +1.62e-03 +1.49e-03 +8.36e-04 +3.95e-04 +2.37e-04 +1.32e-04 +5.59e-05 +1.07e-03 +9.93e-04 +1.63e-03 +6.02e-03 +1.30e-03 +2.86e-03 +2.98e-03 +1.02e-03 +1.31e-03 +1.39e-03 +2.62e-03 +1.69e-03 +2.40e-03 +3.96e-03 +1.26e-03 +1.78e-03 +1.30e-03 +1.06e-03 +1.79e-03 +1.40e-03 +3.18e-03 +1.93e-03 +1.72e-03 +1.78e-03 +8.41e-04 +8.66e-04 +8.00e-04 +1.01e-03 +1.54e-03 +1.30e-03 +1.11e-03 +9.80e-04 +5.02e-04 +3.11e-04 +3.49e-04 +1.24e-03 +1.66e-03 +8.20e-04 +5.29e-04 +3.62e-04 +1.27e-04 +6.03e-05 +6.28e-04 +6.20e-04 +5.01e-04 +4.28e-04 +2.95e-04 +5.38e-05 +6.24e-06 +Upper Bounds +7.61e-01 +7.65e-01 +8.88e-01 +9.87e-01 +1.21e+00 +1.97e+00 +1.13e-02 +7.17e-01 +7.88e-01 +8.74e-01 +9.32e-01 +1.03e+00 +2.41e-01 +1.06e-02 +6.53e-01 +6.86e-01 +8.03e-01 +8.82e-01 +1.02e+00 +5.77e-01 +2.52e-02 +5.23e-01 +6.53e-01 +7.46e-01 +8.58e-01 +9.06e-01 +2.12e-01 +2.66e-02 +3.40e-01 +4.99e-01 +7.82e-01 +8.41e-01 +8.89e-01 +6.94e-02 +2.75e-02 +1.34e-02 +3.33e-02 +1.18e-01 +8.26e-02 +8.54e-02 +5.62e-02 +9.71e-03 +5.62e-04 +2.72e-03 +6.16e-03 +7.54e-03 +7.55e-03 +6.53e-03 +6.05e-03 +7.58e-01 +7.99e-01 +9.04e-01 +1.01e+00 +1.08e+00 +1.29e-01 +1.03e-02 +7.25e-01 +7.95e-01 +9.18e-01 +1.03e+00 +1.07e+00 +4.41e-01 +1.46e-02 +6.80e-01 +7.67e-01 +8.04e-01 +9.61e-01 +1.14e+00 +9.09e-01 +1.35e-02 +5.92e-01 +7.08e-01 +7.60e-01 +9.05e-01 +9.92e-01 +1.43e-01 +3.23e-02 +5.02e-01 +7.06e-01 +9.01e-01 +8.64e-01 +1.02e+00 +5.44e-01 +1.64e-02 +2.68e-02 +7.22e-02 +2.27e-01 +1.08e-01 +3.77e-01 +1.19e-01 +6.91e-03 +2.26e-03 +5.54e-03 +6.95e-03 +6.66e-03 +6.32e-03 +6.08e-03 +4.89e-03 +9.15e-01 +8.87e-01 +9.13e-01 +1.02e+00 +1.13e+00 +3.85e-01 +1.78e-02 +8.13e-01 +9.00e-01 +9.37e-01 +9.28e-01 +1.10e+00 +2.22e-01 +9.99e-03 +8.40e-01 +8.02e-01 +8.17e-01 +8.90e-01 +1.06e+00 +4.93e-01 +1.48e-02 +7.81e-01 +7.93e-01 +8.16e-01 +8.61e-01 +9.05e-01 +3.58e-01 +8.84e-03 +7.91e-01 +9.00e-01 +9.34e-01 +8.31e-01 +8.31e-01 +6.09e-02 +1.33e-02 +1.06e-01 +1.50e-01 +2.43e-01 +7.38e-02 +2.33e-01 +2.69e-02 +5.07e-03 +4.05e-03 +4.28e-03 +7.13e-03 +8.17e-03 +5.26e-03 +5.20e-03 +3.42e-03 +1.01e+00 +1.09e+00 +1.13e+00 +1.18e+00 +1.22e+00 +1.62e-01 +1.34e-01 +1.05e+00 +1.04e+00 +1.03e+00 +1.01e+00 +1.07e+00 +1.21e-01 +1.49e-02 +1.00e+00 +9.68e-01 +9.64e-01 +9.01e-01 +8.62e-01 +6.59e-02 +7.71e-03 +9.19e-01 +9.37e-01 +8.13e-01 +7.43e-01 +7.51e-01 +3.21e-01 +9.55e-03 +9.63e-01 +1.07e+00 +9.67e-01 +7.12e-01 +6.10e-01 +3.60e-02 +4.68e-03 +8.33e-02 +2.66e-01 +2.62e-01 +5.67e-02 +6.16e-02 +2.22e-02 +3.07e-03 +3.63e-03 +3.42e-03 +4.15e-03 +4.88e-03 +3.39e-03 +2.76e-03 +2.65e-03 +1.41e+00 +1.36e+00 +1.42e+00 +1.28e+00 +1.35e+00 +1.12e-01 +1.28e-02 +1.26e+00 +1.24e+00 +1.19e+00 +1.26e+00 +1.18e+00 +7.42e-02 +1.15e-02 +1.09e+00 +1.11e+00 +1.10e+00 +9.79e-01 +9.01e-01 +9.64e-02 +7.63e-03 +1.11e+00 +1.19e+00 +9.95e-01 +8.35e-01 +7.23e-01 +5.53e-02 +4.29e-03 +1.28e+00 +1.16e+00 +9.47e-01 +6.30e-01 +4.16e-01 +4.06e-02 +3.13e-03 +4.27e-01 +3.94e-01 +7.67e-02 +6.12e-02 +3.77e-02 +8.74e-03 +1.70e-03 +5.41e-03 +4.48e-03 +3.06e-03 +2.62e-03 +2.16e-03 +1.25e-03 +1.18e-03 +2.39e+00 +6.90e-01 +2.50e+00 +1.50e-01 +1.56e-01 +5.36e-01 +7.11e-03 +1.50e+00 +2.18e-01 +8.93e-01 +1.20e+00 +3.77e-01 +5.66e-02 +5.84e-03 +1.72e+00 +1.15e+00 +2.77e-01 +2.84e-01 +8.51e-02 +4.01e-02 +4.09e-03 +8.84e-01 +3.28e-01 +3.07e-01 +1.60e-01 +6.60e-02 +3.37e-02 +3.08e-03 +1.74e-01 +8.27e-02 +9.77e-02 +4.36e-02 +4.51e-02 +1.03e-02 +1.71e-03 +3.98e-01 +1.30e-01 +8.25e-02 +3.33e-02 +1.30e-02 +2.87e-03 +7.05e-04 +8.08e-03 +7.43e-03 +4.18e-03 +1.97e-03 +1.18e-03 +6.58e-04 +2.80e-04 +5.36e-03 +4.96e-03 +8.16e-03 +3.01e-02 +6.48e-03 +1.43e-02 +1.49e-02 +5.08e-03 +6.57e-03 +6.95e-03 +1.31e-02 +8.44e-03 +1.20e-02 +1.98e-02 +6.30e-03 +8.90e-03 +6.49e-03 +5.28e-03 +8.94e-03 +7.00e-03 +1.59e-02 +9.64e-03 +8.62e-03 +8.88e-03 +4.21e-03 +4.33e-03 +4.00e-03 +5.04e-03 +7.71e-03 +6.51e-03 +5.56e-03 +4.90e-03 +2.51e-03 +1.55e-03 +1.75e-03 +6.19e-03 +8.31e-03 +4.10e-03 +2.65e-03 +1.81e-03 +6.34e-04 +3.01e-04 +3.14e-03 +3.10e-03 +2.51e-03 +2.14e-03 +1.47e-03 +2.69e-04 +3.12e-05 \ No newline at end of file diff --git a/tests/regression_tests/weightwindows_fw_cadis_local/test.py b/tests/regression_tests/weightwindows_fw_cadis_local/test.py new file mode 100644 index 00000000000..abd3f20986c --- /dev/null +++ b/tests/regression_tests/weightwindows_fw_cadis_local/test.py @@ -0,0 +1,42 @@ +import os + +import openmc +from openmc.examples import random_ray_three_region_cube_with_detectors + +from tests.testing_harness import WeightWindowPyAPITestHarness + + +class MGXSTestHarness(WeightWindowPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_weight_windows_fw_cadis_local(): + model = random_ray_three_region_cube_with_detectors() + + for tally in list(model.tallies): + if tally.name in {"Source Tally", "Absorber Tally", "Cavity Tally"}: + # leave only the tallies of interest + model.tallies.remove(tally) + + ww_mesh = openmc.RegularMesh() + n = 7 + width = 35.0 + ww_mesh.dimension = (n, n, n) + ww_mesh.lower_left = (0.0, 0.0, 0.0) + ww_mesh.upper_right = (width, width, width) + + wwg = openmc.WeightWindowGenerator( + method="fw_cadis", + targets=model.tallies, + mesh=ww_mesh, + max_realizations=model.settings.batches + ) + model.settings.weight_window_generators = wwg + model.settings.random_ray['volume_estimator'] = 'naive' + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat index ceb89e6e34c..a0d84257a8d 100644 --- a/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat @@ -222,11 +222,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true diff --git a/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat index c7691e950c4..62f8478586b 100644 --- a/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat @@ -222,11 +222,13 @@ 500.0 100.0 - - - 0.0 0.0 0.0 30.0 30.0 30.0 - - + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true From aca9dee20c910484aea5296ff9681445652e7a08 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 6 May 2026 16:07:35 +0100 Subject: [PATCH 08/11] Allow Mesh.volumes property for 1D and 2D RegularMesh (#3914) --- openmc/mesh.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 3c3c0a1ac5d..4129913d85c 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -268,10 +268,11 @@ def __repr__(self): return string def _volume_dim_check(self): - if self.n_dimension != 3 or \ - any([d == 0 for d in self.dimension]): - raise RuntimeError(f'Mesh {self.id} is not 3D. ' - 'Volumes cannot be provided.') + if any(d == 0 for d in self.dimension): + raise RuntimeError( + f'Mesh {self.id} has a zero-size dimension. ' + 'Volumes cannot be provided.' + ) @classmethod def from_hdf5(cls, group: h5py.Group): From a4d2d515bf3e6699f0f4f0079b02eed977b4f750 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 7 May 2026 18:50:23 -0500 Subject: [PATCH 09/11] Include mass_attenuation.h5 in package data (#3936) --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index bf40113445b..098487c06ed 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,7 +69,7 @@ include = ['openmc*'] exclude = ['tests*'] [tool.setuptools.package-data] -"openmc.data.dose" = ["**/*.txt"] +"openmc.data.dose" = ["**/*.txt", "*.h5"] "openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"] "openmc.lib" = ["libopenmc.dylib", "libopenmc.so"] From 891a66bc2c1b1afca5737559bba211cfcddf646f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 8 May 2026 20:53:12 -0500 Subject: [PATCH 10/11] Implement DecaySpectrum distribution type and utilize in R2S (#3930) Co-authored-by: Copilot --- docs/source/io_formats/settings.rst | 39 +++- docs/source/pythonapi/stats.rst | 1 + include/openmc/chain.h | 2 + include/openmc/distribution.h | 65 ++++++ openmc/deplete/r2s.py | 246 +++++++++++----------- openmc/stats/univariate.py | 309 ++++++++++++++++++++++++++++ src/chain.cpp | 8 + src/distribution.cpp | 124 +++++++++++ src/finalize.cpp | 2 + src/initialize.cpp | 14 +- src/source.cpp | 30 ++- tests/unit_tests/test_source.py | 69 +++++++ tests/unit_tests/test_stats.py | 135 +++++++++++- 13 files changed, 909 insertions(+), 135 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index d63376e2bdc..f8b064a37d9 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -1058,17 +1058,19 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", - "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces - variates sampled from a uniform distribution over a finite interval. The - "discrete" option produces random variates that can assume a finite number - of values (i.e., a distribution characterized by a probability mass function). - The "tabular" option produces random variates sampled from a tabulated - distribution where the density function is either a histogram or + "tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform" + option produces variates sampled from a uniform distribution over a finite + interval. The "discrete" option produces random variates that can assume a + finite number of values (i.e., a distribution characterized by a probability + mass function). The "tabular" option produces random variates sampled from a + tabulated distribution where the density function is either a histogram or linearly-interpolated between tabulated points. The "watt" option produces random variates is sampled from a Watt fission spectrum (only used for energies). The "maxwell" option produce variates sampled from a Maxwell - fission spectrum (only used for energies). The "mixture" option produces samples - from univariate sub-distributions with given probabilities. + fission spectrum (only used for energies). The "mixture" option produces + samples from univariate sub-distributions with given probabilities. The + "decay_spectrum" option produces photon energies sampled from decay photon + spectra in a depletion chain (only used for energies). *Default*: None @@ -1086,6 +1088,10 @@ variable and whose sub-elements/attributes are as follows: :math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x` points are given first followed by corresponding :math:`p` points. + For a "decay_spectrum" distribution, ``parameters`` gives the atom densities + in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the + same order. + For a "watt" distribution, ``parameters`` should be given as two real numbers :math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c e^{-x/a} \sinh \sqrt{b \, x} dx`. @@ -1115,6 +1121,21 @@ variable and whose sub-elements/attributes are as follows: This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. +:volume: + For a "decay_spectrum" distribution, this attribute specifies the source + region volume in cm\ :sup:`3`. It is used together with atom densities to + determine the absolute photon emission rate. When a source uses a + "decay_spectrum" energy distribution, the source strength is set from this + emission rate. + +:nuclides: + For a "decay_spectrum" distribution, this element specifies a + whitespace-separated list of nuclide names contributing to the decay photon + source. The atom densities for these nuclides are given by the ``parameters`` + element in the same order. Nuclides are resolved against the depletion chain, + and nuclides without decay photon spectra do not contribute to the + distribution. + :bias: This optional element specifies a biased distribution for importance sampling. For continuous distributions, the ``bias`` element should contain another @@ -1300,7 +1321,7 @@ The ```` element specifies the surface flux cosine cutof ```` Element ----------------------------------- -The ```` element specifies the surface flux cosine +The ```` element specifies the surface flux cosine substitution ratio. *Default*: 0.5 diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index 203d4fea4a4..2f2e2835a29 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -22,6 +22,7 @@ Univariate Probability Distributions openmc.stats.Legendre openmc.stats.Mixture openmc.stats.Normal + openmc.stats.DecaySpectrum .. autosummary:: :toctree: generated diff --git a/include/openmc/chain.h b/include/openmc/chain.h index 6f58303585a..e01f2a01c18 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -101,6 +101,8 @@ extern vector> chain_nuclides; void read_chain_file_xml(); +void free_memory_chain(); + } // namespace openmc #endif // OPENMC_CHAIN_H diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index cd81d8801d6..f319dd19df4 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -407,6 +407,71 @@ class Mixture : public Distribution { double integral_; //!< Integral of distribution }; +//============================================================================== +// DecaySpectrum — non-owning mixture of decay photon distributions +//============================================================================== + +//! Energy distribution formed by mixing multiple decay photon spectra. +//! +//! Unlike the general Mixture distribution, this class holds non-owning +//! pointers to the component distributions (which live in +//! data::chain_nuclides). Each component is weighted by the activity +//! (atoms * decay_constant) of the corresponding nuclide. + +class DecaySpectrum : public Distribution { +public: + //============================================================================ + // Types, aliases + + struct Sample { + double energy; + double weight; + int parent_nuclide; + }; + + //============================================================================ + // Constructors + + //! Construct from an XML node containing nuclide names and atom densities. + //! + //! Reads child ```` elements with ``name`` and ``density`` + //! attributes, resolves them against the loaded depletion chain, and + //! constructs the mixed distribution. + explicit DecaySpectrum(pugi::xml_node node); + + //============================================================================ + // Methods + + //! Sample a value from the distribution and return the parent nuclide index + //! \param seed Pseudorandom number seed pointer + //! \return (Sampled energy, sample weight, chain nuclide index) + Sample sample_with_parent(uint64_t* seed) const; + + //! Sample a value from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return (sampled value, sample weight) + std::pair sample(uint64_t* seed) const override; + + double integral() const override; + +protected: + //! Sample a value (unbiased) from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled value + double sample_unbiased(uint64_t* seed) const override; + +private: + //! Initialize decay spectrum sampling data + //! \param nuclide_indices Indices of decay photon emitters in + //! data::chain_nuclides + //! \param atoms Number of atoms for each component. + void init(vector nuclide_indices, const vector& atoms); + + vector nuclide_indices_; //!< Indices of emitting nuclides in the chain + DiscreteIndex di_; //!< Discrete index for component selection + double integral_; //!< Total photon emission rate +}; + } // namespace openmc #endif // OPENMC_DISTRIBUTION_H diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 10d7fdd502a..6dbb3adb2c0 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,5 +1,6 @@ from __future__ import annotations from collections.abc import Sequence +from contextlib import nullcontext import copy from datetime import datetime import json @@ -8,6 +9,7 @@ import numpy as np import openmc from . import IndependentOperator, PredictorIntegrator +from .chain import Chain from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .results import Results from ..checkvalue import PathLike @@ -149,7 +151,7 @@ def run( photon_time_indices: Sequence[int] | None = None, output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, - chain_file: PathLike | None = None, + chain_file: PathLike | Chain | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, @@ -187,9 +189,10 @@ def run( Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method is 'cell-based'. - chain_file : PathLike, optional - Path to the depletion chain XML file to use during activation. If - not provided, the default configured chain file will be used. + chain_file : PathLike or openmc.deplete.Chain, optional + Path to the depletion chain XML file or depletion chain object to + use during activation. If not provided, the default configured + chain file will be used. micro_kwargs : dict, optional Additional keyword arguments passed to :func:`openmc.deplete.get_microxs_and_flux` during the neutron @@ -216,6 +219,8 @@ def run( # consistency (different ranks may have slightly different times) stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S') output_dir = Path(comm.bcast(f'r2s_{stamp}')) + else: + output_dir = Path(output_dir) # Set run_kwargs for the neutron transport step if micro_kwargs is None: @@ -226,22 +231,35 @@ def run( operator_kwargs = {} run_kwargs.setdefault('output', False) micro_kwargs.setdefault('run_kwargs', run_kwargs) - # If a chain file is provided, prefer it for steps 1 and 2 - if chain_file is not None: - micro_kwargs.setdefault('chain_file', chain_file) - operator_kwargs.setdefault('chain_file', chain_file) - self.step1_neutron_transport( - output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs - ) - self.step2_activation( - timesteps, source_rates, timestep_units, output_dir / 'activation', - operator_kwargs=operator_kwargs - ) - self.step3_photon_transport( - photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + # DecaySpectrum distributions are resolved in the C++ solver using + # OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML + # representation alongside the R2S outputs. + if isinstance(chain_file, Chain): + output_dir.mkdir(parents=True, exist_ok=True) + chain_path = output_dir / 'chain.xml' + if comm.rank == 0: + chain_file.export_to_xml(chain_path) + comm.barrier() + else: + chain_path = chain_file + + chain_context = ( + openmc.config.patch('chain_file', chain_path) + if chain_path is not None else nullcontext() ) + with chain_context: + self.step1_neutron_transport( + output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs + ) + self.step2_activation( + timesteps, source_rates, timestep_units, + output_dir / 'activation', operator_kwargs=operator_kwargs + ) + self.step3_photon_transport( + photon_time_indices, bounding_boxes, output_dir / 'photon_transport', + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + ) return output_dir @@ -516,45 +534,30 @@ def step3_photon_transport( if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() - for time_index in time_indices: - # Create decay photon source - if self.method == 'mesh-based': - self.photon_model.settings.source = \ - self.get_decay_photon_source_mesh(time_index) - else: - sources = [] - results = self.results['depletion_results'] - for cell, original_mat in zip(self.domains, self.results['activation_materials']): - # Skip if the cell is not in the photon model or the - # material has changed - if different_photon_model: - if cell.id not in photon_cells or \ + # Determine eligible work items upfront (independent of time index). + if self.method == 'mesh-based': + work_items = self._get_mesh_work_items() + else: + work_items = [] + for cell, original_mat in zip( + self.domains, self.results['activation_materials']): + if different_photon_model: + if cell.id not in photon_cells or \ cell.fill.id != photon_cells[cell.id].fill.id: - continue - - # Get bounding box for the cell - bounding_box = bounding_boxes[cell.id] - - # Get activated material composition - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source source - space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() - strength = energy.integral() if energy is not None else 0.0 - source = openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [cell]} - ) - sources.append(source) - self.photon_model.settings.source = sources + continue + work_items.append((cell, original_mat, bounding_boxes[cell.id])) + + # Ensure photon transport is enabled in settings + self.photon_model.settings.photon_transport = True + for time_index in time_indices: # Convert time_index (which may be negative) to a normal index if time_index < 0: - time_index = len(self.results['depletion_results']) + time_index + time_index += len(self.results['depletion_results']) + + # Build decay photon sources and assign to the photon model + sources = self._create_photon_sources(time_index, work_items) + self.photon_model.settings.source = sources # Run photon transport calculation photon_dir = Path(output_dir) / f'time_{time_index}' @@ -567,58 +570,30 @@ def step3_photon_transport( sp.tallies[tally.id] for tally in self.photon_model.tallies ] - def get_decay_photon_source_mesh( - self, - time_index: int = -1 - ) -> list[openmc.IndependentSource]: - """Create decay photon source for a mesh-based calculation. - - For each mesh element-material combination across all meshes, an - :class:`~openmc.IndependentSource` is created with a - :class:`~openmc.stats.Box` spatial distribution based on the bounding - box of the material within the mesh element. A material constraint is - also applied so that sampled source sites are limited to the correct - region. + def _get_mesh_work_items(self): + """Enumerate mesh-based work items across all meshes. - When the photon transport model is different from the neutron model, the - photon MeshMaterialVolumes is used to determine whether an (element, - material) combination exists in the photon model. - - Parameters - ---------- - time_index : int, optional - Time index for the decay photon source. Default is -1 (last time). + Returns a list of (index_mat, mat_id, bbox) tuples for each eligible + mesh element--material combination, where index_mat is the index into + the activation materials list, mat_id is the material ID, and bbox is + the bounding box for that mesh element--material combination. Returns ------- - list of openmc.IndependentSource - A list of IndependentSource objects for the decay photons, one for - each mesh element-material combination with non-zero source strength. - + list of tuple + Each tuple is (index_mat, mat_id, bbox). """ - mat_dict = self.neutron_model._get_all_materials() - - # List to hold all sources - sources = [] - - # Index in the overall list of activated materials - index_mat = 0 - - # Get various results from previous steps mmv_list = self.results['mesh_material_volumes'] - materials = self.results['activation_materials'] - results = self.results['depletion_results'] photon_mmv_list = self.results.get('mesh_material_volumes_photon') + work_items = [] + index_mat = 0 for mesh_idx, mat_vols in enumerate(mmv_list): photon_mat_vols = photon_mmv_list[mesh_idx] \ if photon_mmv_list is not None else None - # Total number of mesh elements for this mesh n_elements = mat_vols.num_elements - for index_elem in range(n_elements): - # Determine which materials exist in the photon model for this element if photon_mat_vols is not None: photon_materials = { mat_id @@ -626,36 +601,77 @@ def get_decay_photon_source_mesh( if mat_id is not None } - for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): - # Skip void volume + for mat_id, _, bbox in mat_vols.by_element( + index_elem, include_bboxes=True): if mat_id is None: continue - - # Skip if this material doesn't exist in photon model - if photon_mat_vols is not None and mat_id not in photon_materials: + if photon_mat_vols is not None \ + and mat_id not in photon_materials: index_mat += 1 continue - - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source - energy = activated_mat.get_decay_photon_energy() - if energy is not None: - strength = energy.integral() - space = openmc.stats.Box(*bbox) - sources.append(openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [mat_dict[mat_id]]} - )) - - # Increment index of activated material + work_items.append((index_mat, mat_id, bbox)) index_mat += 1 + return work_items + + def _create_photon_sources(self, time_index, work_items): + """Create decay photon sources for a set of regions. + + Builds :class:`openmc.IndependentSource` objects with + :class:`openmc.stats.DecaySpectrum` energy distributions that will be + serialized to XML and resolved against the depletion chain by the C++ + solver. + + Parameters + ---------- + time_index : int + Index into depletion results. + work_items : list of tuple + For mesh-based: list of (index_mat, mat_id, bbox). + For cell-based: list of (cell, original_mat, bbox). + + Returns + ------- + list of openmc.IndependentSource + Photon sources for each activated region. + """ + step_result = self.results['depletion_results'][time_index] + materials = self.results['activation_materials'] + mesh_based = self.method == 'mesh-based' + if mesh_based: + mat_dict = self.neutron_model._get_all_materials() + + sources = [] + for item in work_items: + if mesh_based: + index_mat, domain_id, bbox = item + original_mat = materials[index_mat] + domain = mat_dict[domain_id] + else: + cell, original_mat, bbox = item + domain = cell + + activated_mat = step_result.get_material(str(original_mat.id)) + nuclides = activated_mat.get_nuclide_atom_densities() + if not nuclides: + continue + + # Eliminate nuclides with zero density + nuclides = {nuclide: density for nuclide, density in nuclides.items() + if density > 0} + + energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume) + energy.clip(inplace=True) + if not energy.nuclides: + continue + + sources.append(openmc.IndependentSource( + space=openmc.stats.Box(bbox.lower_left, bbox.upper_right), + energy=energy, + particle='photon', + constraints={'domains': [domain]}, + )) + return sources def load_results(self, path: PathLike): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 1cf9a1ad507..1b18bcb1745 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -2,8 +2,10 @@ from abc import ABC, abstractmethod from collections import defaultdict from collections.abc import Iterable, Sequence +from functools import cache from math import sqrt, pi, exp, log from numbers import Real +from pathlib import Path from warnings import warn import lxml.etree as ET @@ -14,6 +16,7 @@ import openmc.checkvalue as cv from openmc.data import atomic_mass, NEUTRON_MASS +import openmc.data from .._xml import get_elem_list, get_text from ..mixin import EqualityMixin @@ -124,6 +127,8 @@ def from_xml_element(cls, elem): return Legendre.from_xml_element(elem) elif distribution == 'mixture': return Mixture.from_xml_element(elem) + elif distribution == 'decay_spectrum': + return DecaySpectrum.from_xml_element(elem) @abstractmethod def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None): @@ -2196,6 +2201,310 @@ def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: return new_dist +class DecaySpectrum(Univariate): + """Energy distribution from decay photon spectra of a mixture of nuclides. + + This distribution stores nuclide names, their atom densities, and the volume + of the region. When written to XML and read by the C++ solver, the nuclide + names are resolved against the depletion chain to obtain the decay photon + energy spectra and decay constants. The resulting distribution is a mixture + of per-nuclide photon spectra weighted by absolute activity. The volume is + necessary so that the C++ solver can compute the total photon emission rate + in [photons/s], which is used as the source strength. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + nuclides : dict + Dictionary mapping nuclide name (str) to atom density (float) in units + of [atom/b-cm]. + volume : float + Volume of the source region in [cm³]. Used together with atom densities + to compute the absolute photon emission rate. + + Attributes + ---------- + nuclides : dict + Dictionary mapping nuclide name to atom density in [atom/b-cm]. + volume : float + Volume of the source region in [cm³]. + + """ + + def __init__(self, nuclides: dict[str, float], volume: float): + super().__init__(bias=None) + self._dist_cache = None + self._dist_cache_key = None + self.nuclides = nuclides + self.volume = volume + + def __len__(self): + return len(self.nuclides) + + @property + def nuclides(self): + return self._nuclides + + @nuclides.setter + def nuclides(self, nuclides): + cv.check_type('nuclides', nuclides, dict) + for name, density in nuclides.items(): + cv.check_type('nuclide name', name, str) + cv.check_type(f'atom density for {name}', density, Real) + cv.check_greater_than(f'atom density for {name}', density, 0.0, True) + self._nuclides = dict(nuclides) + self._dist_cache = None + self._dist_cache_key = None + + @property + def volume(self): + return self._volume + + @volume.setter + def volume(self, volume): + cv.check_type('volume', volume, Real) + cv.check_greater_than('volume', volume, 0.0) + self._volume = float(volume) + self._dist_cache = None + self._dist_cache_key = None + + @staticmethod + def _chain_file_cache_key(): + """Return a hashable key for the active depletion chain.""" + chain_file = openmc.config.get('chain_file') + if chain_file is None: + return None + + path = Path(chain_file).resolve() + try: + stat = path.stat() + except OSError: + return (path, None, None) + return (path, stat.st_mtime, stat.st_size) + + def to_distribution(self): + """Convert to a concrete distribution using decay chain data. + + Builds a combined photon energy distribution by looking up each nuclide + in the depletion chain via :func:`openmc.data.decay_photon_energy` and + weighting by absolute atom count (``density * 1e24 * volume``). The + result is cached on the object; the cache is invalidated automatically + when :attr:`nuclides` or :attr:`volume` are reassigned. + + Requires ``openmc.config['chain_file']`` to be set. + + Returns + ------- + openmc.stats.Univariate or None + Combined photon energy distribution, or ``None`` if no nuclide in + :attr:`nuclides` has a photon source in the chain. + + """ + chain_key = self._chain_file_cache_key() + if self._dist_cache is not None and self._dist_cache_key == chain_key: + return self._dist_cache + + dists = [] + weights = [] + for name, density in self.nuclides.items(): + dist = openmc.data.decay_photon_energy(name) + if dist is not None: + dists.append(dist) + weights.append(density * 1e24 * self.volume) + + if not dists: + return None + + self._dist_cache = combine_distributions(dists, weights) + self._dist_cache_key = chain_key + return self._dist_cache + + def to_xml_element(self, element_name: str): + """Return XML representation of the decay photon distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : lxml.etree._Element + XML element containing decay photon distribution data + + """ + element = ET.Element(element_name) + element.set("type", "decay_spectrum") + element.set("volume", str(self.volume)) + nuclides = ET.SubElement(element, "nuclides") + nuclides.text = ' '.join(self.nuclides) + parameters = ET.SubElement(element, "parameters") + parameters.text = ' '.join(str(density) for density in self.nuclides.values()) + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element): + """Generate decay photon distribution from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.stats.DecaySpectrum + Decay photon distribution generated from XML element + + """ + volume = float(elem.get('volume')) + names = get_elem_list(elem, 'nuclides', str) + densities = get_elem_list(elem, 'parameters', float) + nuclides = dict(zip(names, densities)) + return cls(nuclides, volume) + + def _sample_unbiased(self, n_samples=1, seed=None): + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum._sample_unbiased requires chain data but none " + "was found. Ensure openmc.config['chain_file'] is set and the " + "chain contains photon sources for the nuclides present." + ) + return dist.sample(n_samples, seed)[0] + + def integral(self): + """Return integral of the distribution + + Returns the total photon emission rate in [photons/s] by delegating to + :meth:`to_distribution`. Returns ``0.0`` when no chain data is + available (e.g., ``openmc.config['chain_file']`` is not set). + + Returns + ------- + float + Total photon emission rate in [photons/s], or ``0.0`` if chain + data is unavailable. + """ + try: + dist = self.to_distribution() + except Exception: + return 0.0 + if dist is None: + return 0.0 + return dist.integral() + + @staticmethod + @cache + def _photon_integral(nuclide: str, chain_key) -> float | None: + """Return the per-atom photon emission integral for a nuclide""" + dist = openmc.data.decay_photon_energy(nuclide) + return dist.integral() if dist is not None else None + + def clip(self, tolerance: float = 1e-9, inplace: bool = False): + """Remove nuclides with negligible contribution to photon emission. + + Nuclides that are stable or have no photon source in the depletion + chain are removed unconditionally. The remaining nuclides are ranked + by their photon emission rate (proportional to + ``atom_density * decay_constant * photon_yield``) and the least + important are discarded until the cumulative discarded fraction of the + total emission rate exceeds *tolerance*. + + Requires ``openmc.config['chain_file']`` to be set. + + Parameters + ---------- + tolerance : float + Maximum fraction of total photon emission rate that may be + discarded. + inplace : bool + Whether to modify the current object in-place or return a new one. + + Returns + ------- + openmc.stats.DecaySpectrum + Distribution with negligible nuclides removed. + + """ + # Compute per-nuclide emission rate; drop non-emitters + emitting_names = [] + emitting_densities = [] + rates = [] + chain_key = self._chain_file_cache_key() + for name, density in self.nuclides.items(): + integral = DecaySpectrum._photon_integral(name, chain_key) + if integral is None: + continue + emitting_names.append(name) + emitting_densities.append(density) + rates.append(density * self.volume * integral) + + if not emitting_names: + new_nuclides = {} + else: + indices = _intensity_clip(rates, tolerance=tolerance) + new_nuclides = { + emitting_names[i]: emitting_densities[i] for i in indices + } + + if inplace: + self._nuclides = new_nuclides + self._dist_cache = None + self._dist_cache_key = None + return self + return type(self)(new_nuclides, self.volume) + + @property + def support(self): + return (0.0, np.inf) + + def evaluate(self, x): + """Evaluate the probability density at a given value. + + Delegates to the combined distribution built from chain data. Raises + ``NotImplementedError`` if the combined distribution is a + :class:`~openmc.stats.Mixture` (which does not support + ``evaluate()``). + + Parameters + ---------- + x : float + Value at which to evaluate the PDF. + + Returns + ------- + float + Probability density at *x*. + """ + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum.evaluate requires chain data. Ensure " + "openmc.config['chain_file'] is set." + ) + return dist.evaluate(x) + + def mean(self): + """Return the mean of the distribution. + + Delegates to the combined distribution built from chain data. + + Returns + ------- + float + Mean photon energy in [eV]. + """ + dist = self.to_distribution() + if dist is None: + raise RuntimeError( + "DecaySpectrum.mean requires chain data. Ensure " + "openmc.config['chain_file'] is set." + ) + return dist.mean() + + def combine_distributions( dists: Sequence[Discrete | Tabular | Mixture], probs: Sequence[float] diff --git a/src/chain.cpp b/src/chain.cpp index e4d0324d3b5..3915c016c29 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -98,6 +98,8 @@ vector> chain_nuclides; void read_chain_file_xml() { + free_memory_chain(); + char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE"); if (!chain_file_path) { return; @@ -120,4 +122,10 @@ void read_chain_file_xml() } } +void free_memory_chain() +{ + data::chain_nuclides.clear(); + data::chain_nuclide_map.clear(); +} + } // namespace openmc diff --git a/src/distribution.cpp b/src/distribution.cpp index c722d0366b6..7f5b498add9 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -7,7 +7,9 @@ #include // for accumulate #include // for runtime_error #include // for string, stod +#include +#include "openmc/chain.h" #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" @@ -15,6 +17,10 @@ #include "openmc/random_lcg.h" #include "openmc/xml_interface.h" +namespace { +std::unordered_set decay_spectrum_missing_chain_nuclides; +} + namespace openmc { //============================================================================== @@ -758,6 +764,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node) dist = UPtrDist {new Tabular(node)}; } else if (type == "mixture") { dist = UPtrDist {new Mixture(node)}; + } else if (type == "decay_spectrum") { + dist = UPtrDist {new DecaySpectrum(node)}; } else if (type == "muir") { openmc::fatal_error( "'muir' distributions are now specified using the openmc.stats.muir() " @@ -768,4 +776,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node) return dist; } +//============================================================================== +// DecaySpectrum implementation +//============================================================================== + +DecaySpectrum::DecaySpectrum(pugi::xml_node node) +{ + // Read the region volume [cm^3] needed for absolute emission rate + if (!check_for_node(node, "volume")) + fatal_error("DecaySpectrum: 'volume' attribute is required."); + double volume = std::stod(get_node_value(node, "volume")); + + // Read nuclide names and atom densities from XML + vector nuclide_indices; + vector atoms; + auto names = get_node_array(node, "nuclides"); + auto densities = get_node_array(node, "parameters"); + if (names.size() != densities.size()) { + fatal_error("DecaySpectrum nuclides and parameters must have the same " + "length."); + } + + for (size_t i = 0; i < names.size(); ++i) { + const auto& name = names[i]; + double density = densities[i]; + + // Look up nuclide in the depletion chain + auto it = data::chain_nuclide_map.find(name); + if (it == data::chain_nuclide_map.end()) { + if (decay_spectrum_missing_chain_nuclides.insert(name).second) { + warning("Nuclide '" + name + + "' appears in a DecaySpectrum source but is not present in " + "the depletion chain; it will be ignored."); + } + continue; + } + + int nuclide_index = it->second; + const auto& chain_nuc = data::chain_nuclides[nuclide_index]; + const Distribution* photon_dist = chain_nuc->photon_energy(); + if (!photon_dist) + continue; + + // Skip non-positive densities and warn if negative + if (density <= 0.0) { + if (density < 0.0) { + warning("Nuclide '" + name + + "' has a negative density in a DecaySpectrum source; it will " + "be ignored."); + } + continue; + } + + // atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3] + double atoms_i = density * 1.0e24 * volume; + + nuclide_indices.push_back(nuclide_index); + atoms.push_back(atoms_i); + } + + init(std::move(nuclide_indices), atoms); +} + +void DecaySpectrum::init( + vector nuclide_indices, const vector& atoms) +{ + if (nuclide_indices.size() != atoms.size()) { + fatal_error("DecaySpectrum nuclide index and atoms arrays must have " + "the same length."); + } + + vector probs; + probs.reserve(nuclide_indices.size()); + for (size_t i = 0; i < nuclide_indices.size(); ++i) { + // Distribution integral is in [photons/s/atom]; multiplying by atoms gives + // the total emission rate [photons/s] for this nuclide. + const auto* dist = + data::chain_nuclides[nuclide_indices[i]]->photon_energy(); + probs.push_back(atoms[i] * dist->integral()); + } + + nuclide_indices_ = std::move(nuclide_indices); + integral_ = std::accumulate(probs.begin(), probs.end(), 0.0); + if (nuclide_indices_.empty() || integral_ <= 0.0) { + fatal_error("DecaySpectrum source did not resolve any nuclides with decay " + "photon spectra and positive atom densities. Ensure " + "OPENMC_CHAIN_FILE is set and matches the nuclides in the " + "source definition."); + } + di_.assign(probs); +} + +DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const +{ + size_t idx = di_.sample(seed); + int parent_nuclide = nuclide_indices_[idx]; + const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy(); + auto [energy, weight] = dist->sample(seed); + return {energy, weight, parent_nuclide}; +} + +std::pair DecaySpectrum::sample(uint64_t* seed) const +{ + auto sample = sample_with_parent(seed); + return {sample.energy, sample.weight}; +} + +double DecaySpectrum::integral() const +{ + return integral_; +} + +double DecaySpectrum::sample_unbiased(uint64_t* seed) const +{ + return sample_with_parent(seed).energy; +} + } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index 98f1125347d..59711dcc31d 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -2,6 +2,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/chain.h" #include "openmc/cmfd_solver.h" #include "openmc/collision_track.h" #include "openmc/constants.h" @@ -44,6 +45,7 @@ void free_memory() free_memory_photon(); free_memory_settings(); free_memory_thermal(); + free_memory_chain(); library_clear(); nuclides_clear(); free_memory_source(); diff --git a/src/initialize.cpp b/src/initialize.cpp index efb462f5c01..991877b2429 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -405,6 +405,10 @@ bool read_model_xml() write_message( fmt::format("Reading model XML file '{}' ...", model_filename), 5); + // Read chain data before settings so DecaySpectrum source distributions can + // resolve nuclides while sources are constructed. + read_chain_file_xml(); + read_settings_xml(settings_root); // If other XML files are present, display warning @@ -420,9 +424,6 @@ bool read_model_xml() } } - // Read data from chain file - read_chain_file_xml(); - // Read materials and cross sections if (!check_for_node(root, "materials")) { fatal_error(fmt::format( @@ -475,14 +476,15 @@ bool read_model_xml() void read_separate_xml_files() { + // Read chain data before settings so DecaySpectrum source distributions can + // resolve nuclides while sources are constructed. + read_chain_file_xml(); + read_settings_xml(); if (settings::run_mode != RunMode::PLOTTING) { read_cross_sections_xml(); } - // Read data from chain file - read_chain_file_xml(); - read_materials_xml(); read_geometry_xml(); diff --git a/src/source.cpp b/src/source.cpp index 524a72bf0e4..bee20d5c305 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -354,6 +354,19 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) if (check_for_node(node, "energy")) { pugi::xml_node node_dist = node.child("energy"); energy_ = distribution_from_xml(node_dist); + + // For decay photon sources, use the absolute photon emission rate in + // [photons/s] as the source strength + if (dynamic_cast(energy_.get())) { + if (strength_ != 1.0) { + warning(fmt::format( + "Source strength of {} is ignored because the source uses a " + "DecaySpectrum energy distribution. The source strength will be " + "set from the DecaySpectrum emission rate.", + strength_)); + } + strength_ = energy_->integral(); + } } else { // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1 energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)}; @@ -414,6 +427,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Check for monoenergetic source above maximum particle energy auto p = particle_.transport_index(); auto energy_ptr = dynamic_cast(energy_.get()); + auto decay_spectrum = dynamic_cast(energy_.get()); if (energy_ptr) { auto energies = tensor::Tensor(energy_ptr->x().data(), energy_ptr->x().size()); @@ -424,10 +438,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const } while (true) { - // Sample energy spectrum - auto [E, E_wgt_temp] = energy_->sample(seed); - site.E = E; - E_wgt = E_wgt_temp; + // Sample energy spectrum. For decay photon sources, also get the parent + // nuclide index to store in the source site for tallying purposes. + if (decay_spectrum) { + auto sample = decay_spectrum->sample_with_parent(seed); + site.E = sample.energy; + E_wgt = sample.weight; + site.parent_nuclide = sample.parent_nuclide; + } else { + auto [E, E_wgt_temp] = energy_->sample(seed); + site.E = E; + E_wgt = E_wgt_temp; + } // Resample if energy falls above maximum particle energy if (site.E < data::energy_max[p] && diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 6eb03f38810..394c09e739d 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -1,5 +1,6 @@ from collections import Counter from math import pi +from pathlib import Path import openmc import openmc.lib @@ -91,6 +92,74 @@ def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model): assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}' +def test_decay_spectrum_parent_nuclide(run_in_tmpdir): + chain_file = Path('chain_decay_spectrum_parent.xml') + chain_file.write_text(""" + + + + 1000000.0 1.0 + + + + + 2000000.0 1.0 + + + +""") + + inner_sphere = openmc.Sphere(r=10.0) + outer_sphere = openmc.Sphere(r=20.0, boundary_type='vacuum') + + shell_mat = openmc.Material() + shell_mat.add_nuclide('H1', 1.0) + shell_mat.set_density('atom/b-cm', 1.0e-12) + + void_cell = openmc.Cell(region=-inner_sphere) + shell_cell = openmc.Cell(fill=shell_mat, region=+inner_sphere & -outer_sphere) + + model = openmc.Model() + model.geometry = openmc.Geometry([void_cell, shell_cell]) + model.materials = [shell_mat] + model.settings.run_mode = 'fixed source' + model.settings.photon_transport = True + model.settings.particles = 1000 + model.settings.batches = 5 + model.settings.source = openmc.IndependentSource( + particle='photon', + space=openmc.stats.Point((0.0, 0.0, 0.0)), + energy=openmc.stats.DecaySpectrum( + {'ParentA': 1.0, 'ParentB': 1.0}, + volume=1.0 + ) + ) + + tally = openmc.Tally() + tally.filters = [ + openmc.CellFilter([void_cell]), + openmc.ParticleFilter(['photon']), + openmc.EnergyFilter([0.0, 1.5e6, 2.5e6]), + openmc.ParentNuclideFilter(['ParentA', 'ParentB']) + ] + tally.scores = ['flux'] + model.tallies = [tally] + + with openmc.config.patch('chain_file', chain_file): + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + tally_out = sp.tallies[tally.id] + mean = tally_out.get_reshaped_data('mean').squeeze() + + assert mean.shape == (2, 2) + assert mean[0, 0] > 0.0 + assert mean[1, 1] > 0.0 + assert mean[0, 1] == 0.0 + assert mean[1, 0] == 0.0 + assert np.count_nonzero(mean) == 2 + + def test_source_file(): filename = 'source.h5' src = openmc.FileSource(path=filename) diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index ca961c8b0f1..2754cf5d09a 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -1,10 +1,11 @@ from math import pi +from pathlib import Path import numpy as np import pytest import openmc import openmc.stats -from openmc.stats.univariate import _INTERPOLATION_SCHEMES +from openmc.stats.univariate import _INTERPOLATION_SCHEMES, DecaySpectrum from scipy.integrate import trapezoid from tests.unit_tests import assert_sample_mean @@ -1013,3 +1014,135 @@ def test_fusion_spectrum_invalid(): # Temperature above 100 keV should raise an error with pytest.raises(ValueError): openmc.stats.fusion_neutron_spectrum(101e3, 'DT') + + +@pytest.fixture(autouse=False) +def decay_spectrum_chain(): + """Set chain_file for the duration of a test and clear the _photon_integral + cache so results from a different chain don't bleed across tests.""" + CHAIN_FILE = (Path(__file__).parents[1] / 'chain_simple.xml').resolve() + DecaySpectrum._photon_integral.cache_clear() + with openmc.config.patch('chain_file', CHAIN_FILE): + yield + DecaySpectrum._photon_integral.cache_clear() + + +def test_decay_spectrum_construction(): + nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4} + d = openmc.stats.DecaySpectrum(nuclides, volume=100.0) + assert d.nuclides == nuclides + assert d.volume == pytest.approx(100.0) + assert len(d) == 2 + + +def test_decay_spectrum_validation(): + # nuclides must be a dict + with pytest.raises(TypeError): + openmc.stats.DecaySpectrum(['I135'], volume=1.0) + + # densities must be > 0 + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': -1.0}, volume=1.0) + + # volume must be > 0 + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=-1.0) + + with pytest.raises(ValueError): + openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=0.0) + + +def test_decay_spectrum_xml_roundtrip(): + nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4} + d = openmc.stats.DecaySpectrum(nuclides, volume=100.0) + + elem = d.to_xml_element('energy') + assert elem.get('type') == 'decay_spectrum' + assert float(elem.get('volume')) == pytest.approx(100.0) + assert elem.findtext('nuclides').split() == list(nuclides) + assert [float(x) for x in elem.findtext('parameters').split()] == pytest.approx( + list(nuclides.values())) + + # Round-trip via DecaySpectrum.from_xml_element + d2 = openmc.stats.DecaySpectrum.from_xml_element(elem) + assert d2.nuclides == nuclides + assert d2.volume == pytest.approx(100.0) + + # Round-trip via the Univariate dispatcher + d3 = openmc.stats.Univariate.from_xml_element(elem) + assert isinstance(d3, openmc.stats.DecaySpectrum) + assert d3 == d + + +def test_decay_spectrum_to_distribution(decay_spectrum_chain): + # Single emitting nuclide -> concrete distribution, not None + d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0) + dist = d.to_distribution() + assert dist is not None + + # Result is cached on second call + dist2 = d.to_distribution() + assert dist2 is dist + + # Nuclide with no photon source -> None + d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + assert d_stable.to_distribution() is None + + # Mixture of emitters -> non-None combined distribution + d_mix = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'Xe135': 5e-4}, volume=10.0 + ) + dist_mix = d_mix.to_distribution() + assert dist_mix is not None + + +def test_decay_spectrum_integral(decay_spectrum_chain): + # For an emitting nuclide, integral should be > 0 + d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0) + assert d.integral() > 0.0 + + # Proportional to density: doubling density doubles integral + d2 = openmc.stats.DecaySpectrum({'I135': 2e-3}, volume=10.0) + assert d2.integral() == pytest.approx(2.0 * d.integral()) + + # Proportional to volume + d3 = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=20.0) + assert d3.integral() == pytest.approx(2.0 * d.integral()) + + # Pure non-emitter -> 0.0 + d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + assert d_stable.integral() == pytest.approx(0.0) + + +def test_decay_spectrum_clip(decay_spectrum_chain): + # Stable / non-emitting nuclides are removed unconditionally + d = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'Xe135': 5e-4, 'Xe136': 1.0, 'Cs135': 1.0}, + volume=10.0, + ) + d_clip = d.clip() + assert 'Xe136' not in d_clip.nuclides + assert 'Cs135' not in d_clip.nuclides + assert 'I135' in d_clip.nuclides + assert 'Xe135' in d_clip.nuclides + # Original is unchanged + assert 'Xe136' in d.nuclides + + # inplace=True modifies and returns the same object + d_same = d.clip(inplace=True) + assert d_same is d + assert 'Xe136' not in d.nuclides + + # A nuclide with negligible emission rate is removed by tolerance clipping. + # U235 has a very small integral (~4e-17 Bq/atom) compared with I135 (~4e-5) + d_tight = openmc.stats.DecaySpectrum( + {'I135': 1e-3, 'U235': 1e-3}, volume=10.0 + ) + d_tight_clip = d_tight.clip(tolerance=1e-9) + assert 'U235' not in d_tight_clip.nuclides + assert 'I135' in d_tight_clip.nuclides + + # All non-emitters -> empty nuclides dict + d_empty = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0) + d_empty.clip(inplace=True) + assert d_empty.nuclides == {} From f76231f92ca5f19948e8e3d62461dbfabe9ed9f8 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Thu, 21 May 2026 13:30:52 -0500 Subject: [PATCH 11/11] fix type error --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 78c2099d28e..0b7cfe4a3bc 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3935,7 +3935,7 @@ void AdaptiveLibMesh::set_mesh_tally_amalgamation( auto elem = *it; auto cluster_elem = elem; - unsigned int cluster_id = + int cluster_id = elem->get_extra_integer(clustering_element_integer_index_); if (cluster_id != -1) {