crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Prevent Adjoint Sources from Trending towards Infinity (#3449)

Browse source
This commit is contained in:
John Tramm 2025-06-18 06:34:56 -05:00 • committed by GitHub
parent 7a1cafa6c1
commit f31130d367
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
2 changed files with 39 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

5
include/openmc/constants.h
View file

@ -63,6 +63,11 @@ constexpr int MAX_SAMPLE {100000};
// source region in the random ray solver
constexpr double MIN_HITS_PER_BATCH {1.5};
// The minimum flux value to be considered non-zero when computing adjoint
// sources. Positive values below this cutoff will be treated as zero, so as to
// prevent extremely large adjoint source terms from being generated.
constexpr double ZERO_FLUX_CUTOFF {1e-22};
// ============================================================================
// MATH AND PHYSICAL CONSTANTS

43
src/random_ray/flat_source_domain.cpp
View file

@ -304,6 +304,13 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
set_flux_to_source(sr, g);
}
}
// Halt if NaN implosion is detected
if (!std::isfinite(source_regions_.scalar_flux_new(sr, g))) {
fatal_error("A source region scalar flux is not finite. "
"This indicates a numerical instability in the "
"simulation. Consider increasing ray density or adjusting "
"the source region mesh.");
}
}
}
@ -1154,9 +1161,9 @@ void FlatSourceDomain::flatten_xs()
m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a);
sigma_t_.push_back(sigma_t);
double nu_Sigma_f =
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);
nu_sigma_f_.push_back(nu_sigma_f);
double sigma_f =
m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
@ -1164,6 +1171,11 @@ void FlatSourceDomain::flatten_xs()
double chi =
m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
if (!std::isfinite(chi)) {
// MGXS interface may return NaN in some cases, such as when material
// is fissionable but has very small sigma_f.
chi = 0.0;
}
chi_.push_back(chi);
for (int g_in = 0; g_in < negroups_; g_in++) {
@ -1191,17 +1203,30 @@ void FlatSourceDomain::flatten_xs()
void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
{
// Set the external source to 1/forward_flux. If the forward flux is negative
// or zero, set the adjoint source to zero, as this is likely a very small
// source region that we don't need to bother trying to vector particles
// towards. Flux negativity in random ray is not related to the flux being
// small in magnitude, but rather due to the source region being physically
// small in volume and thus having a noisy flux estimate.
// 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,
// as this is likely a very small source region that we don't need to bother
// trying to vector particles towards. In the case of flux "being extremely
// close to zero", we define this as being a fixed fraction of the maximum
// forward flux, below which we assume the flux would be physically
// undetectable.
// First, find the maximum forward flux value
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];
if (flux > max_flux) {
max_flux = flux;
}
}
// Then, compute the adjoint source for each source region
#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];
if (flux <= 0.0) {
if (flux <= ZERO_FLUX_CUTOFF * max_flux) {
source_regions_.external_source(sr, g) = 0.0;
} else {
source_regions_.external_source(sr, g) = 1.0 / flux;