crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'develop' into copilot/fix-41682334-bc79-42d2-b03b-308f1a525d17

Browse source
This commit is contained in:
John Tramm 2025-10-02 14:06:31 -05:00
commit 2828a4b7d1
56 changed files with 1089 additions and 473 deletions
Download patch file Download diff file Expand all files Collapse all files

21
docs/source/io_formats/depletion_chain.rst
View file

@ -56,6 +56,27 @@ attributes:
.. _io_chain_reaction:
--------------------
``<source>`` Element
--------------------
The ``<source>`` element represents photon and electron sources associated with
the decay of a nuclide and contains information to construct an
:class:`openmc.stats.Univariate` object that represents this emission as an
energy distribution. This element has the following attributes:
:type:
The type of :class:`openmc.stats.Univariate` source term.
:particle:
The type of particle emitted, e.g., 'photon' or 'electron'
:parameters:
The parameters of the source term, e.g., for a
:class:`openmc.stats.Discrete` source, the energies (in [eV]) at which the
particles are emitted and their relative intensities in [Bq/atom] (in other
words, decay constants).
----------------------
``<reaction>`` Element
----------------------

25
docs/source/usersguide/kinetics.rst
View file

@ -67,6 +67,23 @@ are needed to compute kinetics parameters in OpenMC:
Obtaining kinetics parameters
-----------------------------
The ``Model`` class can be used to automatically generate all IFP tallies using
the Python API with :attr:`openmc.Settings.ifp_n_generation` greater than 0 and
the :meth:`openmc.Model.add_ifp_kinetics_tallies` method::
model = openmc.Model(geometry, settings=settings)
model.add_kinetics_parameters_tallies(num_groups=6) # Add 6 precursor groups
Alternatively, each of the tallies can be manually defined using group-wise or
total :math:`\beta_{\text{eff}}` specified by providing a 6-group
:class:`openmc.DelayedGroupFilter`::
beta_tally = openmc.Tally(name="group-beta-score")
beta_tally.scores = ["ifp-beta-numerator"]
# Add DelayedGroupFilter to enable group-wise tallies
beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))]
Here is an example showing how to declare the three available IFP scores in a
single tally::
@ -95,6 +112,12 @@ for ``ifp-denominator``:
\beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}}
The kinetics parameters can be retrieved directly from a statepoint file using
the :meth:`openmc.StatePoint.ifp_results` method::
with openmc.StatePoint(output_path) as sp:
generation_time, beta_eff = sp.get_kinetics_parameters()
.. only:: html
.. rubric:: References
@ -107,4 +130,4 @@ for ``ifp-denominator``:
of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted
Kinetics Parameters", International Conference on Mathematics and Computational
Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30,
2025 (to be presented).
2025.

5
include/openmc/constants.h
View file

@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5};
// prevent extremely large adjoint source terms from being generated.
constexpr double ZERO_FLUX_CUTOFF {1e-22};
// The minimum macroscopic cross section value considered non-void for the
// random ray solver. Materials with any group with a cross section below this
// value will be converted to pure void.
constexpr double MINIMUM_MACRO_XS {1e-6};
// ============================================================================
// MATH AND PHYSICAL CONSTANTS

7
include/openmc/ifp.h
View file

@ -68,15 +68,14 @@ vector<T> _ifp(const T& value, const vector<T>& data)
//!
//! Add the IFP information in the IFP banks using the same index
//! as the one used to append the fission site to the fission bank.
//! The information stored are the delayed group number and lifetime
//! of the neutron that created the fission event.
//! Multithreading protection is guaranteed by the index returned by the
//! thread_safe_append call in physics.cpp.
//!
//! Needs to be done after the delayed group is found.
//!
//! \param[in] p Particle
//! \param[in] site Fission site
//! \param[in] idx Bank index from the thread_safe_append call in physics.cpp
void ifp(const Particle& p, const SourceSite& site, int64_t idx);
void ifp(const Particle& p, int64_t idx);
//! Resize the IFP banks used in the simulation
void resize_simulation_ifp_banks();

1
include/openmc/particle_data.h
View file

@ -631,6 +631,7 @@ public:
int& event_mt() { return event_mt_; } // MT number of collision
const int& event_mt() const { return event_mt_; }
int& delayed_group() { return delayed_group_; } // delayed group
const int& delayed_group() const { return delayed_group_; }
const int& parent_nuclide() const { return parent_nuclide_; }
int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide

47
include/openmc/random_ray/flat_source_domain.h
View file

@ -27,8 +27,9 @@ public:
//----------------------------------------------------------------------------
// Methods
virtual void update_neutron_source(double k_eff);
double compute_k_eff(double k_eff_old) const;
virtual void update_single_neutron_source(SourceRegionHandle& srh);
virtual void update_all_neutron_sources();
void compute_k_eff();
virtual void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration);
@ -41,7 +42,7 @@ public:
void output_to_vtk() const;
void convert_external_sources();
void count_external_source_regions();
void set_adjoint_sources(const vector<double>& forward_flux);
void set_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;
@ -54,9 +55,8 @@ public:
bool is_target_void);
void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx,
int32_t target_material_id, bool is_target_void);
void prepare_base_source_regions();
SourceRegionHandle get_subdivided_source_region_handle(
int64_t sr, int mesh_bin, Position r, double dist, Direction u);
SourceRegionKey sr_key, Position r, Direction u);
void finalize_discovered_source_regions();
void apply_transport_stabilization();
int64_t n_source_regions() const
@ -67,6 +67,10 @@ public:
{
return source_regions_.n_source_regions() * negroups_;
}
int64_t lookup_base_source_region_idx(const GeometryState& p) const;
SourceRegionKey lookup_source_region_key(const GeometryState& p) const;
int64_t lookup_mesh_bin(int64_t sr, Position r) const;
int lookup_mesh_idx(int64_t sr) const;
//----------------------------------------------------------------------------
// Static Data members
@ -86,6 +90,7 @@ public:
//----------------------------------------------------------------------------
// Public Data members
double k_eff_ {1.0}; // Eigenvalue
bool mapped_all_tallies_ {false}; // If all source regions have been visited
int64_t n_external_source_regions_ {0}; // Total number of source regions with
@ -110,14 +115,6 @@ public:
// The abstract container holding all source region-specific data
SourceRegionContainer source_regions_;
// Base source region container. When source region subdivision via mesh
// is in use, this container holds the original (non-subdivided) material
// filled cell instance source regions. These are useful as they can be
// initialized with external source and mesh domain information ahead of time.
// Then, dynamically discovered source regions can be initialized by cloning
// their base region.
SourceRegionContainer base_source_regions_;
// Parallel hash map holding all source regions discovered during
// a single iteration. This is a threadsafe data structure that is cleaned
// out after each iteration and stored in the "source_regions_" container.
@ -134,8 +131,17 @@ public:
// Map that relates a SourceRegionKey to the external source index. This map
// is used to check if there are any point sources within a subdivided source
// region at the time it is discovered.
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
point_source_map_;
std::unordered_map<SourceRegionKey, vector<int>, SourceRegionKey::HashFunctor>
external_point_source_map_;
// Map that relates a base source region index to the external source index.
// This map is used to check if there are any volumetric sources within a
// subdivided source region at the time it is discovered.
std::unordered_map<int64_t, vector<int>> external_volumetric_source_map_;
// Map that relates a base source region index to a mesh index. This map
// is used to check which subdivision mesh is present in a source region.
std::unordered_map<int64_t, int> mesh_map_;
// If transport corrected MGXS data is being used, there may be negative
// in-group scattering cross sections that can result in instability in MOC
@ -147,12 +153,11 @@ protected:
//----------------------------------------------------------------------------
// Methods
void apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell,
Discrete* discrete, double strength_factor, int target_material_id,
const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(int32_t i_cell,
Discrete* discrete, double strength_factor, int32_t target_material_id);
int src_idx, SourceRegionHandle& srh);
void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx,
int target_material_id, const vector<int32_t>& instances);
void apply_external_source_to_cell_and_children(
int32_t i_cell, int src_idx, int32_t target_material_id);
virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g);
void set_flux_to_source(int64_t sr, int g);
virtual void set_flux_to_old_flux(int64_t sr, int g);

2
include/openmc/random_ray/linear_source_domain.h
View file

@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain {
public:
//----------------------------------------------------------------------------
// Methods
void update_neutron_source(double k_eff) override;
void update_single_neutron_source(SourceRegionHandle& srh) override;
void normalize_scalar_flux_and_volumes(
double total_active_distance_per_iteration) override;

1
include/openmc/random_ray/random_ray.h
View file

@ -48,7 +48,6 @@ public:
static double distance_active_; // Active ray length
static unique_ptr<Source> ray_source_; // Starting source for ray sampling
static RandomRaySourceShape source_shape_; // Flag for linear source
static bool mesh_subdivision_enabled_; // Flag for mesh subdivision
static RandomRaySampleMethod sample_method_; // Flag for sampling method
//----------------------------------------------------------------------------

9
include/openmc/random_ray/random_ray_simulation.h
View file

@ -21,11 +21,7 @@ public:
// Methods
void compute_segment_correction_factors();
void apply_fixed_sources_and_mesh_domains();
void prepare_fixed_sources_adjoint(vector<double>& forward_flux,
SourceRegionContainer& forward_source_regions,
SourceRegionContainer& forward_base_source_regions,
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
forward_source_region_map);
void prepare_fixed_sources_adjoint();
void simulate();
void output_simulation_results() const;
void instability_check(
@ -45,9 +41,6 @@ private:
// Contains all flat source region data
unique_ptr<FlatSourceDomain> domain_;
// Random ray eigenvalue
double k_eff_ {1.0};
// Tracks the average FSR miss rate for analysis and reporting
double avg_miss_rate_ {0.0};

1
include/openmc/random_ray/source_region.h
View file

@ -308,7 +308,6 @@ public:
//----------------------------------------------------------------------------
// Constructors
SourceRegion(int negroups, bool is_linear);
SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
SourceRegion() = default;
//----------------------------------------------------------------------------

2
openmc/data/decay.py
View file

@ -591,7 +591,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None:
openmc.stats.Univariate or None
Distribution of energies in [eV] of photons emitted from decay, or None
if no photon source exists. Note that the probabilities represent
intensities, given as [Bq].
intensities, given as [Bq/atom] (in other words, decay constants).
"""
if not _DECAY_PHOTON_ENERGY:
chain_file = openmc.config.get('chain_file')

13
openmc/mesh.py
View file

@ -287,6 +287,7 @@ class MeshBase(IDManagerMixin, ABC):
model: openmc.Model,
n_samples: int | tuple[int, int, int] = 10_000,
include_void: bool = True,
material_volumes: MeshMaterialVolumes | None = None,
**kwargs
) -> list[openmc.Material]:
"""Generate homogenized materials over each element in a mesh.
@ -305,8 +306,12 @@ class MeshBase(IDManagerMixin, ABC):
the x, y, and z dimensions.
include_void : bool, optional
Whether homogenization should include voids.
material_volumes : MeshMaterialVolumes, optional
Previously computed mesh material volumes to use for homogenization.
If not provided, they will be computed by calling
:meth:`material_volumes`.
**kwargs
Keyword-arguments passed to :meth:`MeshBase.material_volumes`.
Keyword-arguments passed to :meth:`material_volumes`.
Returns
-------
@ -314,7 +319,10 @@ class MeshBase(IDManagerMixin, ABC):
Homogenized material in each mesh element
"""
vols = self.material_volumes(model, n_samples, **kwargs)
if material_volumes is None:
vols = self.material_volumes(model, n_samples, **kwargs)
else:
vols = material_volumes
mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
# Create homogenized material for each element
@ -424,7 +432,6 @@ class MeshBase(IDManagerMixin, ABC):
# Restore original tallies
model.tallies = original_tallies
return volumes

41
openmc/model/model.py
View file

@ -1,7 +1,7 @@
from __future__ import annotations
from collections.abc import Iterable, Sequence
import copy
from functools import lru_cache
from functools import cache
from pathlib import Path
import math
from numbers import Integral, Real
@ -160,7 +160,7 @@ class Model:
return False
@property
@lru_cache(maxsize=None)
@cache
def _materials_by_id(self) -> dict:
"""Dictionary mapping material ID --> material"""
if self.materials:
@ -170,14 +170,14 @@ class Model:
return {mat.id: mat for mat in mats}
@property
@lru_cache(maxsize=None)
@cache
def _cells_by_id(self) -> dict:
"""Dictionary mapping cell ID --> cell"""
cells = self.geometry.get_all_cells()
return {cell.id: cell for cell in cells.values()}
@property
@lru_cache(maxsize=None)
@cache
def _cells_by_name(self) -> dict[int, openmc.Cell]:
# Get the names maps, but since names are not unique, store a set for
# each name key. In this way when the user requests a change by a name,
@ -190,7 +190,7 @@ class Model:
return result
@property
@lru_cache(maxsize=None)
@cache
def _materials_by_name(self) -> dict[int, openmc.Material]:
if self.materials is None:
mats = self.geometry.get_all_materials().values()
@ -203,6 +203,37 @@ class Model:
result[mat.name].add(mat)
return result
def add_kinetics_parameters_tallies(self, num_groups: int | None = None):
"""Add tallies for calculating kinetics parameters using the IFP method.
This method adds tallies to the model for calculating two kinetics
parameters, the generation time and the effective delayed neutron
fraction (beta effective). After a model is run, these parameters can be
determined through the :meth:`openmc.StatePoint.ifp_results` method.
Parameters
----------
num_groups : int, optional
Number of precursor groups to filter the delayed neutron fraction.
If None, only the total effective delayed neutron fraction is
tallied.
"""
if not any('ifp-time-numerator' in t.scores for t in self.tallies):
gen_time_tally = openmc.Tally(name='IFP time numerator')
gen_time_tally.scores = ['ifp-time-numerator']
self.tallies.append(gen_time_tally)
if not any('ifp-beta-numerator' in t.scores for t in self.tallies):
beta_tally = openmc.Tally(name='IFP beta numerator')
beta_tally.scores = ['ifp-beta-numerator']
if num_groups is not None:
beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))]
self.tallies.append(beta_tally)
if not any('ifp-denominator' in t.scores for t in self.tallies):
denom_tally = openmc.Tally(name='IFP denominator')
denom_tally.scores = ['ifp-denominator']
self.tallies.append(denom_tally)
@classmethod
def from_xml(
cls,

67
openmc/statepoint.py
View file

@ -1,4 +1,5 @@
from datetime import datetime
from collections import namedtuple
import glob
import re
import os
@ -8,6 +9,7 @@ import h5py
import numpy as np
from pathlib import Path
from uncertainties import ufloat
from uncertainties.unumpy import uarray
import openmc
import openmc.checkvalue as cv
@ -15,6 +17,9 @@ import openmc.checkvalue as cv
_VERSION_STATEPOINT = 18
KineticsParameters = namedtuple("KineticsParameters", ["generation_time", "beta_effective"])
class StatePoint:
"""State information on a simulation at a certain point in time (at the end
of a given batch). Statepoints can be used to analyze tally results as well
@ -531,7 +536,7 @@ class StatePoint:
def get_tally(self, scores=[], filters=[], nuclides=[],
name=None, id=None, estimator=None, exact_filters=False,
exact_nuclides=False, exact_scores=False,
multiply_density=None, derivative=None):
multiply_density=None, derivative=None, filter_type=None):
"""Finds and returns a Tally object with certain properties.
This routine searches the list of Tallies and returns the first Tally
@ -575,6 +580,9 @@ class StatePoint:
to the same value as this parameter.
derivative : openmc.TallyDerivative, optional
TallyDerivative object to match.
filter_type : type, optional
If not None, the Tally must have at least one Filter that is an
instance of this type. For example `openmc.MeshFilter`.
Returns
-------
@ -648,6 +656,10 @@ class StatePoint:
if not contains_filters:
continue
if filter_type is not None:
if not any(isinstance(f, filter_type) for f in test_tally.filters):
continue
# Determine if Tally has the queried Nuclide(s)
if nuclides:
if not all(nuclide in test_tally.nuclides for nuclide in nuclides):
@ -710,3 +722,56 @@ class StatePoint:
tally_filter.paths = cell.paths
self._summary = summary
def get_kinetics_parameters(self) -> KineticsParameters:
"""Get kinetics parameters from IFP tallies.
This method searches the tallies in the statepoint for the tallies
required to compute kinetics parameters using the Iterated Fission
Probability (IFP) method.
Returns
-------
KineticsParameters
A named tuple containing the generation time and effective delayed
neutron fraction. If the necessary tallies for one or both
parameters are not found, that parameter is returned as None.
"""
denom_tally = None
gen_time_tally = None
beta_tally = None
for tally in self.tallies.values():
if 'ifp-denominator' in tally.scores:
denom_tally = self.get_tally(scores=['ifp-denominator'])
if 'ifp-time-numerator' in tally.scores:
gen_time_tally = self.get_tally(scores=['ifp-time-numerator'])
if 'ifp-beta-numerator' in tally.scores:
beta_tally = self.get_tally(scores=['ifp-beta-numerator'])
if denom_tally is None:
return KineticsParameters(None, None)
def get_ufloat(tally, score):
return uarray(tally.get_values(scores=[score]),
tally.get_values(scores=[score], value='std_dev'))
denom_values = get_ufloat(denom_tally, 'ifp-denominator')
if gen_time_tally is None:
generation_time = None
else:
gen_time_values = get_ufloat(gen_time_tally, 'ifp-time-numerator')
gen_time_values /= denom_values*self.keff
generation_time = gen_time_values.flatten()[0]
if beta_tally is None:
beta_effective = None
else:
beta_values = get_ufloat(beta_tally, 'ifp-beta-numerator')
beta_values /= denom_values
beta_effective = beta_values.flatten()
if beta_effective.size == 1:
beta_effective = beta_effective[0]
return KineticsParameters(generation_time, beta_effective)

4
src/ifp.cpp
View file

@ -28,13 +28,13 @@ bool is_generation_time_or_both()
return false;
}
void ifp(const Particle& p, const SourceSite& site, int64_t idx)
void ifp(const Particle& p, int64_t idx)
{
if (is_beta_effective_or_both()) {
const auto& delayed_groups =
simulation::ifp_source_delayed_group_bank[p.current_work() - 1];
simulation::ifp_fission_delayed_group_bank[idx] =
_ifp(site.delayed_group, delayed_groups);
_ifp(p.delayed_group(), delayed_groups);
}
if (is_generation_time_or_both()) {
const auto& lifetimes =

1
src/particle.cpp
View file

@ -144,6 +144,7 @@ void Particle::from_source(const SourceSite* src)
time() = src->time;
time_last() = src->time;
parent_nuclide() = src->parent_nuclide;
delayed_group() = src->delayed_group;
// Convert signed surface ID to signed index
if (src->surf_id != SURFACE_NONE) {

9
src/physics.cpp
View file

@ -246,18 +246,15 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
}
// Iterated Fission Probability (IFP) method
if (settings::ifp_on) {
ifp(p, site, idx);
ifp(p, idx);
}
} else {
p.secondary_bank().push_back(site);
}
// Set the delayed group on the particle as well
p.delayed_group() = site.delayed_group;
// Increment the number of neutrons born delayed
if (p.delayed_group() > 0) {
nu_d[p.delayed_group() - 1]++;
if (site.delayed_group > 0) {
nu_d[site.delayed_group - 1]++;
}
// Write fission particles to nuBank

373
src/random_ray/flat_source_domain.cpp
View file

@ -53,24 +53,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
// Initialize source regions.
bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
source_regions_ = SourceRegionContainer(negroups_, is_linear);
source_regions_.assign(
base_source_regions, SourceRegion(negroups_, is_linear));
// Initialize materials
int64_t source_region_id = 0;
for (int i = 0; i < model::cells.size(); i++) {
Cell& cell = *model::cells[i];
if (cell.type_ == Fill::MATERIAL) {
for (int j = 0; j < cell.n_instances(); j++) {
source_regions_.material(source_region_id++) = cell.material(j);
}
}
}
// Sanity check
if (source_region_id != base_source_regions) {
fatal_error("Unexpected number of source regions");
}
// Initialize tally volumes
if (volume_normalized_flux_tallies_) {
@ -118,34 +100,24 @@ void FlatSourceDomain::accumulate_iteration_flux()
}
}
// Compute new estimate of scattering + fission sources in each source region
// based on the flux estimate from the previous iteration.
void FlatSourceDomain::update_neutron_source(double k_eff)
void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
{
simulation::time_update_src.start();
double inverse_k_eff = 1.0 / k_eff;
// Reset all source regions to zero (important for void regions)
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 0.0;
// Reset all source regions to zero (important for void regions)
for (int g = 0; g < negroups_; g++) {
srh.source(g) = 0.0;
}
// Add scattering + fission source
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
int material = srh.material();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material * negroups_ + g_out];
double scatter_source = 0.0;
double fission_source = 0.0;
for (int g_in = 0; g_in < negroups_; g_in++) {
double scalar_flux = source_regions_.scalar_flux_old(sr, g_in);
double scalar_flux = srh.scalar_flux_old(g_in);
double sigma_s =
sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];
double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];
@ -154,18 +126,30 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
scatter_source += sigma_s * scalar_flux;
fission_source += nu_sigma_f * scalar_flux * chi;
}
source_regions_.source(sr, g_out) =
srh.source(g_out) =
(scatter_source + fission_source * inverse_k_eff) / sigma_t;
}
}
// Add external source if in fixed source mode
if (settings::run_mode == RunMode::FIXED_SOURCE) {
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
for (int g = 0; g < negroups_; g++) {
srh.source(g) += srh.external_source(g);
}
}
}
// Compute new estimate of scattering + fission sources in each source region
// based on the flux estimate from the previous iteration.
void FlatSourceDomain::update_all_neutron_sources()
{
simulation::time_update_src.start();
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
update_single_neutron_source(srh);
}
simulation::time_update_src.stop();
}
@ -320,7 +304,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
// Generates new estimate of k_eff based on the differences between this
// iteration's estimate of the scalar flux and the last iteration's estimate.
double FlatSourceDomain::compute_k_eff(double k_eff_old) const
void FlatSourceDomain::compute_k_eff()
{
double fission_rate_old = 0;
double fission_rate_new = 0;
@ -365,7 +349,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
p[sr] = sr_fission_source_new;
}
double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old);
double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old);
double H = 0.0;
// defining an inverse sum for better performance
@ -385,7 +369,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
// Adds entropy value to shared entropy vector in openmc namespace.
simulation::entropy.push_back(H);
return k_eff_new;
k_eff_ = k_eff_new;
}
// This function is responsible for generating a mapping between random
@ -652,7 +636,6 @@ void FlatSourceDomain::random_ray_tally()
"random ray mode.");
break;
}
// Apply score to the appropriate tally bin
Tally& tally {*model::tallies[task.tally_idx]};
#pragma omp atomic
@ -726,21 +709,21 @@ void FlatSourceDomain::output_to_vtk() const
print_plot();
// Outer loop over plots
for (int p = 0; p < model::plots.size(); p++) {
for (int plt = 0; plt < model::plots.size(); plt++) {
// Get handle to OpenMC plot object and extract params
Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[p].get());
Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[plt].get());
// Random ray plots only support voxel plots
if (!openmc_plot) {
warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting "
"is allowed in random ray mode.",
p));
plt));
continue;
} else if (openmc_plot->type_ != Plot::PlotType::voxel) {
warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting "
"is allowed in random ray mode.",
p));
plt));
continue;
}
@ -794,23 +777,11 @@ void FlatSourceDomain::output_to_vtk() const
continue;
}
int i_cell = p.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + p.cell_instance();
if (RandomRay::mesh_subdivision_enabled_) {
int mesh_idx = base_source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
mesh_bin = model::meshes[mesh_idx]->get_bin(p.r());
}
SourceRegionKey sr_key {sr, mesh_bin};
auto it = source_region_map_.find(sr_key);
if (it != source_region_map_.end()) {
sr = it->second;
} else {
sr = -1;
}
SourceRegionKey sr_key = lookup_source_region_key(p);
int64_t sr = -1;
auto it = source_region_map_.find(sr_key);
if (it != source_region_map_.end()) {
sr = it->second;
}
voxel_indices[z * Ny * Nx + y * Nx + x] = sr;
@ -967,13 +938,17 @@ void FlatSourceDomain::output_to_vtk() const
}
void FlatSourceDomain::apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, SourceRegionHandle& srh)
int src_idx, SourceRegionHandle& srh)
{
srh.external_source_present() = 1;
auto s = model::external_sources[src_idx].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto discrete = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();
const auto& discrete_energies = discrete->x();
const auto& discrete_probs = discrete->prob();
srh.external_source_present() = 1;
for (int i = 0; i < discrete_energies.size(); i++) {
int g = data::mg.get_group_index(discrete_energies[i]);
srh.external_source(g) += discrete_probs[i] * strength_factor;
@ -981,8 +956,7 @@ void FlatSourceDomain::apply_external_source_to_source_region(
}
void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
Discrete* discrete, double strength_factor, int target_material_id,
const vector<int32_t>& instances)
int src_idx, int target_material_id, const vector<int32_t>& instances)
{
Cell& cell = *model::cells[i_cell];
@ -1000,16 +974,13 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell,
if (target_material_id == C_NONE ||
cell_material_id == target_material_id) {
int64_t source_region = source_region_offsets_[i_cell] + j;
SourceRegionHandle srh =
source_regions_.get_source_region_handle(source_region);
apply_external_source_to_source_region(discrete, strength_factor, srh);
external_volumetric_source_map_[source_region].push_back(src_idx);
}
}
}
void FlatSourceDomain::apply_external_source_to_cell_and_children(
int32_t i_cell, Discrete* discrete, double strength_factor,
int32_t target_material_id)
int32_t i_cell, int src_idx, int32_t target_material_id)
{
Cell& cell = *model::cells[i_cell];
@ -1017,14 +988,14 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children(
vector<int> instances(cell.n_instances());
std::iota(instances.begin(), instances.end(), 0);
apply_external_source_to_cell_instances(
i_cell, discrete, strength_factor, target_material_id, instances);
i_cell, src_idx, target_material_id, instances);
} else if (target_material_id == C_NONE) {
std::unordered_map<int32_t, vector<int32_t>> cell_instance_list =
cell.get_contained_cells(0, nullptr);
for (const auto& pair : cell_instance_list) {
int32_t i_child_cell = pair.first;
apply_external_source_to_cell_instances(i_child_cell, discrete,
strength_factor, target_material_id, pair.second);
apply_external_source_to_cell_instances(
i_child_cell, src_idx, target_material_id, pair.second);
}
}
}
@ -1070,36 +1041,17 @@ void FlatSourceDomain::convert_external_sources()
"point source at {}",
sp->r()));
}
int i_cell = gs.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
SourceRegionKey key = lookup_source_region_key(gs);
if (RandomRay::mesh_subdivision_enabled_) {
// If mesh subdivision is enabled, we need to determine which subdivided
// mesh bin the point source coordinate is in as well
int mesh_idx = source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r());
}
// With the source region and mesh bin known, we can use the
// accompanying SourceRegionKey as a key into a map that stores the
// corresponding external source index for the point source. Notably, we
// do not actually apply the external source to any source regions here,
// as if mesh subdivision is enabled, they haven't actually been
// discovered & initilized yet. When discovered, they will read from the
// point_source_map to determine if there are any point source terms
// that should be applied.
SourceRegionKey key {sr, mesh_bin};
point_source_map_[key] = es;
} else {
// If we are not using mesh subdivision, we can apply the external
// source directly to the source region as we do for volumetric domain
// constraint sources.
SourceRegionHandle srh = source_regions_.get_source_region_handle(sr);
apply_external_source_to_source_region(energy, strength_factor, srh);
}
// With the source region and mesh bin known, we can use the
// accompanying SourceRegionKey as a key into a map that stores the
// corresponding external source index for the point source. Notably, we
// do not actually apply the external source to any source regions here,
// as if mesh subdivision is enabled, they haven't actually been
// discovered & initilized yet. When discovered, they will read from the
// external_source_map to determine if there are any external source
// terms that should be applied.
external_point_source_map_[key].push_back(es);
} else {
// If not a point source, then use the volumetric domain constraints to
@ -1107,42 +1059,25 @@ void FlatSourceDomain::convert_external_sources()
if (is->domain_type() == Source::DomainType::MATERIAL) {
for (int32_t material_id : domain_ids) {
for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, material_id);
apply_external_source_to_cell_and_children(i_cell, es, material_id);
}
}
} else if (is->domain_type() == Source::DomainType::CELL) {
for (int32_t cell_id : domain_ids) {
int32_t i_cell = model::cell_map[cell_id];
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, C_NONE);
apply_external_source_to_cell_and_children(i_cell, es, C_NONE);
}
} else if (is->domain_type() == Source::DomainType::UNIVERSE) {
for (int32_t universe_id : domain_ids) {
int32_t i_universe = model::universe_map[universe_id];
Universe& universe = *model::universes[i_universe];
for (int32_t i_cell : universe.cells_) {
apply_external_source_to_cell_and_children(
i_cell, energy, strength_factor, C_NONE);
apply_external_source_to_cell_and_children(i_cell, es, C_NONE);
}
}
}
}
} // End loop over external sources
// Divide the fixed source term by sigma t (to save time when applying each
// iteration)
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
source_regions_.external_source(sr, g) /= sigma_t;
}
}
}
void FlatSourceDomain::flux_swap()
@ -1159,13 +1094,23 @@ void FlatSourceDomain::flatten_xs()
const int a = 0;
n_materials_ = data::mg.macro_xs_.size();
for (auto& m : data::mg.macro_xs_) {
for (int i = 0; i < n_materials_; i++) {
auto& m = data::mg.macro_xs_[i];
for (int g_out = 0; g_out < negroups_; g_out++) {
if (m.exists_in_model) {
double sigma_t =
m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
sigma_t_.push_back(sigma_t);
if (sigma_t < MINIMUM_MACRO_XS) {
Material* mat = model::materials[i].get();
warning(fmt::format(
"Material \"{}\" (id: {}) has a group {} total cross section "
"({:.3e}) below the minimum threshold "
"({:.3e}). Material will be treated as pure void.",
mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS));
}
double nu_sigma_f =
m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a);
nu_sigma_f_.push_back(nu_sigma_f);
@ -1206,7 +1151,7 @@ void FlatSourceDomain::flatten_xs()
}
}
void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
void FlatSourceDomain::set_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,
@ -1220,7 +1165,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
double max_flux = 0.0;
#pragma omp parallel for reduction(max : max_flux)
for (int64_t se = 0; se < n_source_elements(); se++) {
double flux = forward_flux[se];
double flux = source_regions_.scalar_flux_final(se);
if (flux > max_flux) {
max_flux = flux;
}
@ -1230,7 +1175,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
for (int g = 0; g < negroups_; g++) {
double flux = forward_flux[sr * negroups_ + g];
double flux = source_regions_.scalar_flux_final(sr, g);
if (flux <= ZERO_FLUX_CUTOFF * max_flux) {
source_regions_.external_source(sr, g) = 0.0;
} else {
@ -1239,6 +1184,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
if (flux > 0.0) {
source_regions_.external_source_present(sr) = 1;
}
source_regions_.scalar_flux_final(sr, g) = 0.0;
}
}
@ -1265,7 +1211,6 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
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
@ -1326,13 +1271,14 @@ void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell,
if ((target_material_id == C_NONE && !is_target_void) ||
cell_material_id == target_material_id) {
int64_t sr = source_region_offsets_[i_cell] + j;
if (source_regions_.mesh(sr) != C_NONE) {
// print out the source region that is broken:
// Check if the key is already present in the mesh_map_
if (mesh_map_.find(sr) != mesh_map_.end()) {
fatal_error(fmt::format("Source region {} already has mesh idx {} "
"applied, but trying to apply mesh idx {}",
sr, source_regions_.mesh(sr), mesh_idx));
sr, mesh_map_[sr], mesh_idx));
}
source_regions_.mesh(sr) = mesh_idx;
// If the SR has not already been assigned, then we can write to it
mesh_map_[sr] = mesh_idx;
}
}
}
@ -1402,18 +1348,9 @@ void FlatSourceDomain::apply_meshes()
}
}
void FlatSourceDomain::prepare_base_source_regions()
{
std::swap(source_regions_, base_source_regions_);
source_regions_.negroups() = base_source_regions_.negroups();
source_regions_.is_linear() = base_source_regions_.is_linear();
}
SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
int64_t sr, int mesh_bin, Position r, double dist, Direction u)
SourceRegionKey sr_key, Position r, Direction u)
{
SourceRegionKey sr_key {sr, mesh_bin};
// Case 1: Check if the source region key is already present in the permanent
// map. This is the most common condition, as any source region visited in a
// previous power iteration will already be present in the permanent map. If
@ -1475,9 +1412,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
gs.r() = r + TINY_BIT * u;
gs.u() = {1.0, 0.0, 0.0};
exhaustive_find_cell(gs);
int gs_i_cell = gs.lowest_coord().cell();
int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance();
if (sr_found != sr) {
int64_t sr_found = lookup_base_source_region_idx(gs);
if (sr_found != sr_key.base_source_region_id) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1485,9 +1421,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
}
// Sanity check on mesh bin
int mesh_idx = base_source_regions_.mesh(sr);
int mesh_idx = lookup_mesh_idx(sr_key.base_source_region_id);
if (mesh_idx == C_NONE) {
if (mesh_bin != 0) {
if (sr_key.mesh_bin != 0) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1496,7 +1432,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
int bin_found = mesh->get_bin(r + TINY_BIT * u);
if (bin_found != mesh_bin) {
if (bin_found != sr_key.mesh_bin) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
@ -1508,26 +1444,60 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
// condition only occurs the first time the source region is discovered
// (typically in the first power iteration). In this case, we need to handle
// creation of the new source region and its storage into the parallel map.
// The new source region is created by copying the base source region, so as
// to inherit material, external source, and some flux properties etc. We
// also pass the base source region id to allow the new source region to
// know which base source region it is derived from.
SourceRegion* sr_ptr = discovered_source_regions_.emplace(
sr_key, {base_source_regions_.get_source_region_handle(sr), sr});
discovered_source_regions_.unlock(sr_key);
// Additionally, we need to determine the source region's material, initialize
// the starting scalar flux guess, and apply any known external sources.
// Call the basic constructor for the source region and store in the parallel
// map.
bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
SourceRegion* sr_ptr =
discovered_source_regions_.emplace(sr_key, {negroups_, is_linear});
SourceRegionHandle handle {*sr_ptr};
// Check if the new source region contains a point source and apply it if so
auto it2 = point_source_map_.find(sr_key);
if (it2 != point_source_map_.end()) {
int es = it2->second;
auto s = model::external_sources[es].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto energy = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();
apply_external_source_to_source_region(energy, strength_factor, handle);
int material = handle.material();
if (material != MATERIAL_VOID) {
// Determine the material
int gs_i_cell = gs.lowest_coord().cell();
Cell& cell = *model::cells[gs_i_cell];
int material = cell.material(gs.cell_instance());
// If material total XS is extremely low, just set it to void to avoid
// problems with 1/Sigma_t
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
if (sigma_t < MINIMUM_MACRO_XS) {
material = MATERIAL_VOID;
break;
}
}
handle.material() = material;
// Store the mesh index (if any) assigned to this source region
handle.mesh() = mesh_idx;
if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Determine if there are any volumetric sources, and apply them.
// Volumetric sources are specifc only to the base SR idx.
auto it_vol =
external_volumetric_source_map_.find(sr_key.base_source_region_id);
if (it_vol != external_volumetric_source_map_.end()) {
const vector<int>& vol_sources = it_vol->second;
for (int src_idx : vol_sources) {
apply_external_source_to_source_region(src_idx, handle);
}
}
// Determine if there are any point sources, and apply them.
// Point sources are specific to the source region key.
auto it_point = external_point_source_map_.find(sr_key);
if (it_point != external_point_source_map_.end()) {
const vector<int>& point_sources = it_point->second;
for (int src_idx : point_sources) {
apply_external_source_to_source_region(src_idx, handle);
}
}
// Divide external source term by sigma_t
if (material != C_NONE) {
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
handle.external_source(g) /= sigma_t;
@ -1535,6 +1505,21 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
}
}
// Compute the combined source term
update_single_neutron_source(handle);
// Unlock the parallel map. Note: we may be tempted to release
// this lock earlier, and then just use the source region's lock to protect
// the flux/source initialization stages above. However, the rest of the code
// only protects updates to the new flux and volume fields, and assumes that
// the source is constant for the duration of transport. Thus, using just the
// source region's lock by itself would result in other threads potentially
// reading from the source before it is computed, as they won't use the lock
// when only reading from the SR's source. It would be expensive to protect
// those operations, whereas generating the SR is only done once, so we just
// hold the map's bucket lock until the source region is fully initialized.
discovered_source_regions_.unlock(sr_key);
return handle;
}
@ -1620,4 +1605,52 @@ void FlatSourceDomain::apply_transport_stabilization()
}
}
// Determines the base source region index (i.e., a material filled cell
// instance) that corresponds to a particular location in the geometry. Requires
// that the "gs" object passed in has already been initialized and has called
// find_cell etc.
int64_t FlatSourceDomain::lookup_base_source_region_idx(
const GeometryState& gs) const
{
int i_cell = gs.lowest_coord().cell();
int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
return sr;
}
// Determines the index of the mesh (if any) that has been applied
// to a particular base source region index.
int FlatSourceDomain::lookup_mesh_idx(int64_t sr) const
{
int mesh_idx = C_NONE;
auto mesh_it = mesh_map_.find(sr);
if (mesh_it != mesh_map_.end()) {
mesh_idx = mesh_it->second;
}
return mesh_idx;
}
// Determines the source region key that corresponds to a particular location in
// the geometry. This takes into account both the base source region index as
// well as the mesh bin if a mesh is applied to this source region for
// subdivision.
SourceRegionKey FlatSourceDomain::lookup_source_region_key(
const GeometryState& gs) const
{
int64_t sr = lookup_base_source_region_idx(gs);
int64_t mesh_bin = lookup_mesh_bin(sr, gs.r());
return SourceRegionKey {sr, mesh_bin};
}
// Determines the mesh bin that corresponds to a particular base source region
// index and position.
int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const
{
int mesh_idx = lookup_mesh_idx(sr);
int mesh_bin = 0;
if (mesh_idx != C_NONE) {
mesh_bin = model::meshes[mesh_idx]->get_bin(r);
}
return mesh_bin;
}
} // namespace openmc

47
src/random_ray/linear_source_domain.cpp
View file

@ -34,25 +34,18 @@ void LinearSourceDomain::batch_reset()
}
}
void LinearSourceDomain::update_neutron_source(double k_eff)
void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
{
simulation::time_update_src.start();
double inverse_k_eff = 1.0 / k_eff;
// Reset all source regions to zero (important for void regions)
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 0.0;
// Reset all source regions to zero (important for void regions)
for (int g = 0; g < negroups_; g++) {
srh.source(g) = 0.0;
}
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
}
MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();
// Add scattering + fission source
int material = srh.material();
if (material != MATERIAL_VOID) {
double inverse_k_eff = 1.0 / k_eff_;
MomentMatrix invM = srh.mom_matrix().inverse();
for (int g_out = 0; g_out < negroups_; g_out++) {
double sigma_t = sigma_t_[material * negroups_ + g_out];
@ -64,8 +57,8 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
for (int g_in = 0; g_in < negroups_; g_in++) {
// Handles for the flat and linear components of the flux
double flux_flat = source_regions_.scalar_flux_old(sr, g_in);
MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in);
double flux_flat = srh.scalar_flux_old(g_in);
MomentArray flux_linear = srh.flux_moments_old(g_in);
// Handles for cross sections
double sigma_s =
@ -81,7 +74,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
}
// Compute the flat source term
source_regions_.source(sr, g_out) =
srh.source(g_out) =
(scatter_flat + fission_flat * inverse_k_eff) / sigma_t;
// Compute the linear source terms. In the first 10 iterations when the
@ -91,25 +84,21 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
// very small/noisy or have poorly developed spatial moments, so we zero
// the source gradients (effectively making this a flat source region
// temporarily), so as to improve stability.
if (simulation::current_batch > 10 &&
source_regions_.source(sr, g_out) >= 0.0) {
source_regions_.source_gradients(sr, g_out) =
if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) {
srh.source_gradients(g_out) =
invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
} else {
source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0};
srh.source_gradients(g_out) = {0.0, 0.0, 0.0};
}
}
}
// Add external source if in fixed source mode
if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Add external source to flat source term if in fixed source mode
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
for (int g = 0; g < negroups_; g++) {
srh.source(g) += srh.external_source(g);
}
}
simulation::time_update_src.stop();
}
void LinearSourceDomain::normalize_scalar_flux_and_volumes(

121
src/random_ray/random_ray.cpp
View file

@ -237,7 +237,6 @@ double RandomRay::distance_inactive_;
double RandomRay::distance_active_;
unique_ptr<Source> RandomRay::ray_source_;
RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT};
bool RandomRay::mesh_subdivision_enabled_ {false};
RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG};
RandomRay::RandomRay()
@ -336,71 +335,60 @@ void RandomRay::event_advance_ray()
void RandomRay::attenuate_flux(double distance, bool is_active, double offset)
{
// Determine source region index etc.
int i_cell = lowest_coord().cell();
// The base source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
// Lookup base source region index
int64_t sr = domain_->lookup_base_source_region_idx(*this);
// Perform ray tracing across mesh
if (mesh_subdivision_enabled_) {
// Determine the mesh index for the base source region, if any
int mesh_idx = domain_->base_source_regions_.mesh(sr);
// Determine the mesh index for the base source region, if any
int mesh_idx = domain_->lookup_mesh_idx(sr);
if (mesh_idx == C_NONE) {
// If there's no mesh being applied to this cell, then
// we just attenuate the flux as normal, and set
// the mesh bin to 0
attenuate_flux_inner(distance, is_active, sr, 0, r());
} else {
// If there is a mesh being applied to this cell, then
// we loop over all the bin crossings and attenuate
// separately.
Mesh* mesh = model::meshes[mesh_idx].get();
// We adjust the start and end positions of the ray slightly
// to accomodate for floating point precision issues that tend
// to occur at mesh boundaries that overlap with geometry lattice
// boundaries.
Position start = r() + (offset + TINY_BIT) * u();
Position end = start + (distance - 2.0 * TINY_BIT) * u();
double reduced_distance = (end - start).norm();
// Ray trace through the mesh and record bins and lengths
mesh_bins_.resize(0);
mesh_fractional_lengths_.resize(0);
mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_);
// Loop over all mesh bins and attenuate flux
for (int b = 0; b < mesh_bins_.size(); b++) {
double physical_length = reduced_distance * mesh_fractional_lengths_[b];
attenuate_flux_inner(
physical_length, is_active, sr, mesh_bins_[b], start);
start += physical_length * u();
}
}
if (mesh_idx == C_NONE) {
// If there's no mesh being applied to this cell, then
// we just attenuate the flux as normal, and set
// the mesh bin to 0
attenuate_flux_inner(distance, is_active, sr, 0, r());
} else {
attenuate_flux_inner(distance, is_active, sr, C_NONE, r());
// If there is a mesh being applied to this cell, then
// we loop over all the bin crossings and attenuate
// separately.
Mesh* mesh = model::meshes[mesh_idx].get();
// We adjust the start and end positions of the ray slightly
// to accomodate for floating point precision issues that tend
// to occur at mesh boundaries that overlap with geometry lattice
// boundaries.
Position start = r() + (offset + TINY_BIT) * u();
Position end = start + (distance - 2.0 * TINY_BIT) * u();
double reduced_distance = (end - start).norm();
// Ray trace through the mesh and record bins and lengths
mesh_bins_.resize(0);
mesh_fractional_lengths_.resize(0);
mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_);
// Loop over all mesh bins and attenuate flux
for (int b = 0; b < mesh_bins_.size(); b++) {
double physical_length = reduced_distance * mesh_fractional_lengths_[b];
attenuate_flux_inner(
physical_length, is_active, sr, mesh_bins_[b], start);
start += physical_length * u();
}
}
}
void RandomRay::attenuate_flux_inner(
double distance, bool is_active, int64_t sr, int mesh_bin, Position r)
{
SourceRegionKey sr_key {sr, mesh_bin};
SourceRegionHandle srh;
if (mesh_subdivision_enabled_) {
srh = domain_->get_subdivided_source_region_handle(
sr, mesh_bin, r, distance, u());
if (srh.is_numerical_fp_artifact_) {
return;
}
} else {
srh = domain_->source_regions_.get_source_region_handle(sr);
srh = domain_->get_subdivided_source_region_handle(sr_key, r, u());
if (srh.is_numerical_fp_artifact_) {
return;
}
switch (source_shape_) {
case RandomRaySourceShape::FLAT:
if (this->material() == MATERIAL_VOID) {
if (srh.material() == MATERIAL_VOID) {
attenuate_flux_flat_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_flat_source(srh, distance, is_active, r);
@ -408,7 +396,7 @@ void RandomRay::attenuate_flux_inner(
break;
case RandomRaySourceShape::LINEAR:
case RandomRaySourceShape::LINEAR_XY:
if (this->material() == MATERIAL_VOID) {
if (srh.material() == MATERIAL_VOID) {
attenuate_flux_linear_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_linear_source(srh, distance, is_active, r);
@ -439,7 +427,7 @@ void RandomRay::attenuate_flux_flat_source(
n_event()++;
// Get material
int material = this->material();
int material = srh.material();
// MOC incoming flux attenuation + source contribution/attenuation equation
for (int g = 0; g < negroups_; g++) {
@ -490,7 +478,7 @@ void RandomRay::attenuate_flux_flat_source_void(
// The number of geometric intersections is counted for reporting purposes
n_event()++;
int material = this->material();
int material = srh.material();
// If ray is in the active phase (not in dead zone), make contributions to
// source region bookkeeping
@ -537,7 +525,7 @@ void RandomRay::attenuate_flux_linear_source(
// The number of geometric intersections is counted for reporting purposes
n_event()++;
int material = this->material();
int material = srh.material();
Position& centroid = srh.centroid();
Position midpoint = r + u() * (distance / 2.0);
@ -810,27 +798,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
cell_born() = lowest_coord().cell();
}
SourceRegionKey sr_key = domain_->lookup_source_region_key(*this);
SourceRegionHandle srh =
domain_->get_subdivided_source_region_handle(sr_key, r(), u());
// Initialize ray's starting angular flux to starting location's isotropic
// source
int i_cell = lowest_coord().cell();
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
SourceRegionHandle srh;
if (mesh_subdivision_enabled_) {
int mesh_idx = domain_->base_source_regions_.mesh(sr);
int mesh_bin;
if (mesh_idx == C_NONE) {
mesh_bin = 0;
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
mesh_bin = mesh->get_bin(r());
}
srh =
domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u());
} else {
srh = domain_->source_regions_.get_source_region_handle(sr);
}
if (!srh.is_numerical_fp_artifact_) {
for (int g = 0; g < negroups_; g++) {
angular_flux_[g] = srh.source(g);

190
src/random_ray/random_ray_simulation.cpp
View file

@ -47,97 +47,82 @@ void openmc_run_random_ray()
if (mpi::master)
validate_random_ray_inputs();
// Declare forward flux so that it can be saved for later adjoint simulation
vector<double> forward_flux;
SourceRegionContainer forward_source_regions;
SourceRegionContainer forward_base_source_regions;
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
forward_source_region_map;
// Initialize Random Ray Simulation Object
RandomRaySimulation sim;
{
// Initialize Random Ray Simulation Object
RandomRaySimulation sim;
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
// Begin main simulation timer
simulation::time_total.start();
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
sim.simulate();
// Execute random ray simulation
sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// End main simulation timer
simulation::time_total.stop();
// Normalize and save the final forward flux
sim.domain()->serialize_final_fluxes(forward_flux);
double source_normalization_factor =
sim.domain()->compute_fixed_source_normalization_factor() /
(settings::n_batches - settings::n_inactive);
// Normalize and save the final forward flux
double source_normalization_factor =
sim.domain()->compute_fixed_source_normalization_factor() /
(settings::n_batches - settings::n_inactive);
#pragma omp parallel for
for (uint64_t i = 0; i < forward_flux.size(); i++) {
forward_flux[i] *= source_normalization_factor;
}
forward_source_regions = sim.domain()->source_regions_;
forward_source_region_map = sim.domain()->source_region_map_;
forward_base_source_regions = sim.domain()->base_source_regions_;
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) {
sim.domain()->source_regions_.scalar_flux_final(se) *=
source_normalization_factor;
}
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
//////////////////////////////////////////////////////////
// Run adjoint simulation (if enabled)
//////////////////////////////////////////////////////////
if (adjoint_needed) {
reset_timers();
// Configure the domain for adjoint simulation
FlatSourceDomain::adjoint_ = true;
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
// Initialize OpenMC general data structures
openmc_simulation_init();
// Initialize Random Ray Simulation Object
RandomRaySimulation adjoint_sim;
// Initialize adjoint fixed sources, if present
adjoint_sim.prepare_fixed_sources_adjoint(forward_flux,
forward_source_regions, forward_base_source_regions,
forward_source_region_map);
// Transpose scattering matrix
adjoint_sim.domain()->transpose_scattering_matrix();
// Swap nu_sigma_f and chi
adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_);
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
adjoint_sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
adjoint_sim.output_simulation_results();
if (!adjoint_needed) {
return;
}
reset_timers();
// Configure the domain for adjoint simulation
FlatSourceDomain::adjoint_ = true;
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
// Initialize OpenMC general data structures
openmc_simulation_init();
sim.domain()->k_eff_ = 1.0;
// Initialize adjoint fixed sources, if present
sim.prepare_fixed_sources_adjoint();
// Transpose scattering matrix
sim.domain()->transpose_scattering_matrix();
// Swap nu_sigma_f and chi
sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_);
// Begin main simulation timer
simulation::time_total.start();
// Execute random ray simulation
sim.simulate();
// End main simulation timer
simulation::time_total.stop();
// Finalize OpenMC
openmc_simulation_finalize();
// Output all simulation results
sim.output_simulation_results();
}
// Enforces restrictions on inputs in random ray mode. While there are
@ -348,7 +333,6 @@ void validate_random_ray_inputs()
// when generating weight windows with FW-CADIS and an overlaid mesh.
///////////////////////////////////////////////////////////////////
if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR &&
RandomRay::mesh_subdivision_enabled_ &&
variance_reduction::weight_windows.size() > 0) {
warning(
"Linear sources may result in negative fluxes in small source regions "
@ -366,7 +350,6 @@ void openmc_reset_random_ray()
FlatSourceDomain::mesh_domain_map_.clear();
RandomRay::ray_source_.reset();
RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
RandomRay::mesh_subdivision_enabled_ = false;
RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
}
@ -412,20 +395,11 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
}
}
void RandomRaySimulation::prepare_fixed_sources_adjoint(
vector<double>& forward_flux, SourceRegionContainer& forward_source_regions,
SourceRegionContainer& forward_base_source_regions,
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
forward_source_region_map)
void RandomRaySimulation::prepare_fixed_sources_adjoint()
{
domain_->source_regions_.adjoint_reset();
if (settings::run_mode == RunMode::FIXED_SOURCE) {
if (RandomRay::mesh_subdivision_enabled_) {
domain_->source_regions_ = forward_source_regions;
domain_->source_region_map_ = forward_source_region_map;
domain_->base_source_regions_ = forward_base_source_regions;
domain_->source_regions_.adjoint_reset();
}
domain_->set_adjoint_sources(forward_flux);
domain_->set_adjoint_sources();
}
}
@ -445,22 +419,18 @@ void RandomRaySimulation::simulate()
simulation::total_weight = 1.0;
// Update source term (scattering + fission)
domain_->update_neutron_source(k_eff_);
domain_->update_all_neutron_sources();
// Reset scalar fluxes, iteration volume tallies, and region hit flags to
// zero
// Reset scalar fluxes, iteration volume tallies, and region hit flags
// to zero
domain_->batch_reset();
// At the beginning of the simulation, if mesh subvivision is in use, we
// At the beginning of the simulation, if mesh subdivision is in use, we
// need to swap the main source region container into the base container,
// as the main source region container will be used to hold the true
// subdivided source regions. The base container will therefore only
// contain the external source region information, the mesh indices,
// material properties, and initial guess values for the flux/source.
if (RandomRay::mesh_subdivision_enabled_ &&
simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) {
domain_->prepare_base_source_regions();
}
// Start timer for transport
simulation::time_transport.start();
@ -476,11 +446,9 @@ void RandomRaySimulation::simulate()
simulation::time_transport.stop();
// If using mesh subdivision, add any newly discovered source regions
// to the main source region container.
if (RandomRay::mesh_subdivision_enabled_) {
domain_->finalize_discovered_source_regions();
}
// Add any newly discovered source regions to the main source region
// container.
domain_->finalize_discovered_source_regions();
// Normalize scalar flux and update volumes
domain_->normalize_scalar_flux_and_volumes(
@ -494,10 +462,10 @@ void RandomRaySimulation::simulate()
if (settings::run_mode == RunMode::EIGENVALUE) {
// Compute random ray k-eff
k_eff_ = domain_->compute_k_eff(k_eff_);
domain_->compute_k_eff();
// Store random ray k-eff into OpenMC's native k-eff variable
global_tally_tracklength = k_eff_;
global_tally_tracklength = domain_->k_eff_;
}
// Execute all tallying tasks, if this is an active batch
@ -507,12 +475,6 @@ void RandomRaySimulation::simulate()
// estimate
domain_->accumulate_iteration_flux();
// Generate mapping between source regions and tallies
if (!domain_->mapped_all_tallies_ &&
!RandomRay::mesh_subdivision_enabled_) {
domain_->convert_source_regions_to_tallies(0);
}
// Use above mapping to contribute FSR flux data to appropriate
// tallies
domain_->random_ray_tally();
@ -522,7 +484,7 @@ void RandomRaySimulation::simulate()
domain_->flux_swap();
// Check for any obvious insabilities/nans/infs
instability_check(n_hits, k_eff_, avg_miss_rate_);
instability_check(n_hits, domain_->k_eff_, avg_miss_rate_);
} // End MPI master work
// Finalize the current batch
@ -571,7 +533,7 @@ void RandomRaySimulation::instability_check(
}
if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) {
fatal_error("Instability detected");
fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff));
}
}
}

28
src/random_ray/source_region.cpp
View file

@ -48,7 +48,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear)
}
scalar_flux_new_.assign(negroups, 0.0);
source_.resize(negroups);
source_.assign(negroups, 0.0);
scalar_flux_final_.assign(negroups, 0.0);
tally_task_.resize(negroups);
@ -60,25 +60,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear)
}
}
SourceRegion::SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr)
: SourceRegion(handle.negroups_, handle.is_linear_)
{
material_ = handle.material();
mesh_ = handle.mesh();
parent_sr_ = parent_sr;
for (int g = 0; g < scalar_flux_new_.size(); g++) {
scalar_flux_old_[g] = handle.scalar_flux_old(g);
source_[g] = handle.source(g);
}
if (settings::run_mode == RunMode::FIXED_SOURCE) {
external_source_present_ = handle.external_source_present();
for (int g = 0; g < scalar_flux_new_.size(); g++) {
external_source_[g] = handle.external_source(g);
}
}
}
//==============================================================================
// SourceRegionContainer implementation
//==============================================================================
@ -259,9 +240,12 @@ void SourceRegionContainer::adjoint_reset()
MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(),
MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0);
if (settings::run_mode == RunMode::FIXED_SOURCE) {
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0);
} else {
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0);
}
std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0);
std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0);
std::fill(source_.begin(), source_.end(), 0.0f);
std::fill(external_source_.begin(), external_source_.end(), 0.0f);
std::fill(source_gradients_.begin(), source_gradients_.end(),

1
src/settings.cpp
View file

@ -346,7 +346,6 @@ void get_run_parameters(pugi::xml_node node_base)
}
FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back(
type, domain_id);
RandomRay::mesh_subdivision_enabled_ = true;
}
}
}

3
src/tallies/tally.cpp
View file

@ -560,7 +560,8 @@ void Tally::set_scores(const vector<std::string>& scores)
// Make sure a delayed group filter wasn't used with an incompatible
// score.
if (delayedgroup_filter_ != C_NONE) {
if (score_str != "delayed-nu-fission" && score_str != "decay-rate")
if (score_str != "delayed-nu-fission" && score_str != "decay-rate" &&
score_str != "ifp-beta-numerator")
fatal_error("Cannot tally " + score_str + "with a delayedgroup filter");
}

9
src/tallies/tally_scoring.cpp
View file

@ -964,6 +964,15 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
if (delayed_groups.size() == settings::ifp_n_generation) {
if (delayed_groups[0] > 0) {
score = p.wgt_last();
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
const DelayedGroupFilter& filt {
*dynamic_cast<DelayedGroupFilter*>(
model::tally_filters[i_dg_filt].get())};
score_fission_delayed_dg(i_tally, delayed_groups[0] - 1,
score, score_index, p.filter_matches());
continue;
}
}
}
}

0
tests/regression_tests/ifp/groupwise/__init__.py Normal file
View file

43
tests/regression_tests/ifp/groupwise/inputs_true.dat Normal file
View file

@ -0,0 +1,43 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1" name="core" depletable="true">
<density value="16.0" units="g/cm3"/>
<nuclide name="U235" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" region="-1" universe="1"/>
<surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 10.0"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>20</batches>
<inactive>5</inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-10.0 -10.0 -10.0 10.0 10.0 10.0</parameters>
</space>
<constraints>
<fissionable>true</fissionable>
</constraints>
</source>
<ifp_n_generation>5</ifp_n_generation>
</settings>
<tallies>
<filter id="1" type="delayedgroup">
<bins>1 2 3 4 5 6</bins>
</filter>
<tally id="1" name="IFP time numerator">
<scores>ifp-time-numerator</scores>
</tally>
<tally id="2" name="IFP beta numerator">
<filters>1</filters>
<scores>ifp-beta-numerator</scores>
</tally>
<tally id="3" name="IFP denominator">
<scores>ifp-denominator</scores>
</tally>
</tallies>
</model>

21
tests/regression_tests/ifp/groupwise/results_true.dat Normal file
View file

@ -0,0 +1,21 @@
k-combined:
1.006559E+00 5.389391E-03
tally 1:
9.109384E-08
5.667165E-16
tally 2:
3.000000E-03
9.000000E-06
0.000000E+00
0.000000E+00
2.100000E-02
1.370000E-04
2.800000E-02
2.220000E-04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
tally 3:
1.489000E+01
1.480036E+01

40
tests/regression_tests/ifp/groupwise/test.py Normal file
View file

@ -0,0 +1,40 @@
"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted
kinetics parameters using dedicated tallies."""
import openmc
import pytest
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def ifp_model():
# Material
material = openmc.Material(name="core")
material.add_nuclide("U235", 1.0)
material.set_density('g/cm3', 16.0)
# Geometry
radius = 10.0
sphere = openmc.Sphere(r=radius, boundary_type="vacuum")
cell = openmc.Cell(region=-sphere, fill=material)
geometry = openmc.Geometry([cell])
# Settings
settings = openmc.Settings()
settings.particles = 1000
settings.batches = 20
settings.inactive = 5
settings.ifp_n_generation = 5
model = openmc.Model(settings=settings, geometry=geometry)
space = openmc.stats.Box(*cell.bounding_box)
model.settings.source = openmc.IndependentSource(
space=space, constraints={'fissionable': True})
model.add_kinetics_parameters_tallies(num_groups=6)
return model
def test_iterated_fission_probability(ifp_model):
harness = PyAPITestHarness("statepoint.20.h5", model=ifp_model)
harness.main()

0
tests/regression_tests/ifp/total/__init__.py Normal file
View file

0
tests/regression_tests/ifp/inputs_true.dat → tests/regression_tests/ifp/total/inputs_true.dat
View file

4
tests/regression_tests/ifp/results_true.dat → tests/regression_tests/ifp/total/results_true.dat
View file

@ -3,7 +3,7 @@ k-combined:
tally 1:
9.109384E-08
5.667165E-16
6.500000E-02
6.710000E-04
5.200000E-02
5.420000E-04
1.489000E+01
1.480036E+01

1
tests/regression_tests/ifp/test.py → tests/regression_tests/ifp/total/test.py
View file

@ -6,7 +6,6 @@ import pytest
from tests.testing_harness import PyAPITestHarness
@pytest.fixture()
def ifp_model():
model = openmc.Model()

0
tests/regression_tests/random_ray_low_density/__init__.py Normal file
View file

244
tests/regression_tests/random_ray_low_density/inputs_true.dat Normal file
View file

@ -0,0 +1,244 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="source">
<density value="1.0" units="macro"/>
<macroscopic name="source"/>
</material>
<material id="2" name="void">
<density value="1.0" units="macro"/>
<macroscopic name="void"/>
</material>
<material id="3" name="absorber">
<density value="1.0" units="macro"/>
<macroscopic name="absorber"/>
</material>
</materials>
<geometry>
<cell id="1" name="infinite source region" material="1" universe="1"/>
<cell id="2" name="infinite void region" material="2" universe="2"/>
<cell id="3" name="infinite absorber region" material="3" universe="3"/>
<cell id="4" fill="4" universe="5"/>
<cell id="5" name="full domain" fill="5" region="1 -2 3 -4 5 -6" universe="6"/>
<lattice id="4">
<pitch>2.5 2.5 2.5</pitch>
<dimension>12 12 12</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<universes>
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 </universes>
</lattice>
<surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
<surface id="2" type="x-plane" boundary="vacuum" coeffs="30.0"/>
<surface id="3" type="y-plane" boundary="reflective" coeffs="0.0"/>
<surface id="4" type="y-plane" boundary="vacuum" coeffs="30.0"/>
<surface id="5" type="z-plane" boundary="reflective" coeffs="0.0"/>
<surface id="6" type="z-plane" boundary="vacuum" coeffs="30.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>90</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="3.14" particle="neutron">
<energy type="discrete">
<parameters>100.0 1.0</parameters>
</energy>
<constraints>
<domain_type>universe</domain_type>
<domain_ids>1</domain_ids>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>500.0</distance_active>
<distance_inactive>100.0</distance_inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
</random_ray>
</settings>
<tallies>
<filter id="3" type="material">
<bins>1</bins>
</filter>
<filter id="2" type="material">
<bins>2</bins>
</filter>
<filter id="1" type="material">
<bins>3</bins>
</filter>
<tally id="3" name="Source Tally">
<filters>3</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="2" name="Void Tally">
<filters>2</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="1" name="Absorber Tally">
<filters>1</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>
</model>

9
tests/regression_tests/random_ray_low_density/results_true.dat Normal file
View file

@ -0,0 +1,9 @@
tally 1:
5.973607E-01
7.155477E-02
tally 2:
3.206216E-02
2.063375E-04
tally 3:
2.096415E-03
8.804963E-07

60
tests/regression_tests/random_ray_low_density/test.py Normal file
View file

@ -0,0 +1,60 @@
import os
import numpy as np
import openmc
from openmc.examples import random_ray_three_region_cube
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_low_density():
model = random_ray_three_region_cube()
# Rebuild the MGXS library to have a material with very
# low macroscopic cross sections
ebins = [1e-5, 20.0e6]
groups = openmc.mgxs.EnergyGroups(group_edges=ebins)
void_sigma_a = 4.0e-6
void_sigma_s = 3.0e-4
void_mat_data = openmc.XSdata('void', groups)
void_mat_data.order = 0
void_mat_data.set_total([void_sigma_a + void_sigma_s])
void_mat_data.set_absorption([void_sigma_a])
void_mat_data.set_scatter_matrix(
np.rollaxis(np.array([[[void_sigma_s]]]), 0, 3))
absorber_sigma_a = 0.75
absorber_sigma_s = 0.25
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.0000001
source_sigma_a = void_sigma_a * multiplier
source_sigma_s = void_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, void_mat_data, absorber_mat_data])
mg_cross_sections_file.export_to_hdf5()
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()

8
tests/regression_tests/random_ray_point_source_locator/results_true.dat
View file

@ -1,9 +1,9 @@
tally 1:
2.633900E+00
2.948207E+00
2.633923E+00
2.948228E+00
tally 2:
1.440463E-01
3.294032E-03
1.440456E-01
3.293984E-03
tally 3:
9.425207E-03
1.089748E-05

BIN
tests/regression_tests/surface_source_write/case-03/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-04/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-05/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-07/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-08/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-09/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-10/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-12/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-13/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-14/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-d07/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-d08/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-e02/surface_source_true.h5
View file

Binary file not shown.

BIN
tests/regression_tests/surface_source_write/case-e03/surface_source_true.h5
View file

Binary file not shown.

39
tests/unit_tests/test_ifp.py
View file

@ -47,3 +47,42 @@ def test_exceptions(options, error, run_in_tmpdir, geometry):
tallies = openmc.Tallies([tally])
model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies)
model.run()
@pytest.mark.parametrize(
"num_groups, use_auto_tallies",
[
(None, True),
(None, False),
(6, True),
(6, False),
],
)
def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_tallies):
# Create basic model
model = openmc.Model(geometry=geometry)
model.settings.particles = 1000
model.settings.batches = 20
model.settings.inactive = 5
model.settings.ifp_n_generation = 5
# Add IFP tallies either via the convenience method or manually
if use_auto_tallies:
model.add_kinetics_parameters_tallies(num_groups=num_groups)
else:
for score in ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"]:
tally = openmc.Tally()
tally.scores = [score]
if score == "ifp-beta-numerator" and num_groups is not None:
tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))]
model.tallies.append(tally)
# Run and get kinetics parameters
sp_file = model.run()
with openmc.StatePoint(sp_file) as sp:
params = sp.get_kinetics_parameters()
assert isinstance(params, openmc.KineticsParameters)
assert params.generation_time is not None
assert params.beta_effective is not None
if num_groups is not None:
assert len(params.beta_effective) == num_groups

65
tests/unit_tests/test_statepoint.py Normal file
View file

@ -0,0 +1,65 @@
import openmc
def test_get_tally_filter_type(run_in_tmpdir):
"""Test various ways of retrieving tallies from a StatePoint object."""
mat = openmc.Material()
mat.add_nuclide("H1", 1.0)
mat.set_density("g/cm3", 10.0)
sphere = openmc.Sphere(r=10.0, boundary_type="vacuum")
cell = openmc.Cell(fill=mat, region=-sphere)
geometry = openmc.Geometry([cell])
settings = openmc.Settings()
settings.particles = 10
settings.batches = 2
settings.run_mode = "fixed source"
reg_mesh = openmc.RegularMesh().from_domain(cell)
tally1 = openmc.Tally(tally_id=1)
mesh_filter = openmc.MeshFilter(reg_mesh)
tally1.filters = [mesh_filter]
tally1.scores = ["flux"]
tally2 = openmc.Tally(tally_id=2, name="heating tally")
cell_filter = openmc.CellFilter(cell)
tally2.filters = [cell_filter]
tally2.scores = ["heating"]
tallies = openmc.Tallies([tally1, tally2])
model = openmc.Model(
geometry=geometry, materials=[mat], settings=settings, tallies=tallies
)
sp_filename = model.run()
sp = openmc.StatePoint(sp_filename)
tally_found = sp.get_tally(filter_type=openmc.MeshFilter)
assert tally_found.id == 1
tally_found = sp.get_tally(filter_type=openmc.CellFilter)
assert tally_found.id == 2
tally_found = sp.get_tally(filters=[mesh_filter])
assert tally_found.id == 1
tally_found = sp.get_tally(filters=[cell_filter])
assert tally_found.id == 2
tally_found = sp.get_tally(scores=["heating"])
assert tally_found.id == 2
tally_found = sp.get_tally(name="heating tally")
assert tally_found.id == 2
tally_found = sp.get_tally(name=None)
assert tally_found.id == 1
tally_found = sp.get_tally(id=1)
assert tally_found.id == 1
tally_found = sp.get_tally(id=2)
assert tally_found.id == 2