Random Ray Point Source Locator (#3360)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
John Tramm 2025-04-02 00:40:07 -05:00 committed by GitHub
parent b67771a3d0
commit 07f5334616
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 435 additions and 33 deletions

View file

@ -1052,7 +1052,8 @@ random ray and Monte Carlo, however.
regions. Thus, in the OpenMC implementation of random ray, particle sources
are restricted to being volumetric and isotropic, although different energy
spectrums are supported. Fixed sources can be applied to specific materials,
cells, or universes.
cells, or universes. Point sources are "smeared" to fill the volume of the
source region that contains the point source coordinate.
- **Inactive batches:** In Monte Carlo, use of a fixed source implies that all
batches are active batches, as there is no longer a need to develop a fission

View file

@ -919,9 +919,12 @@ Monte Carlo solver.
Currently, all of the following conditions must be met for the particle source
to be valid in random ray mode:
- One or more domain ids must be specified that indicate which cells, universes,
or materials the source applies to. This implicitly limits the source type to
being volumetric. This is specified via the ``domains`` constraint placed on the
- Either a point source must be used, or a domain constraint must be specified
that indicates which cells, universes, or materials the source applies to. In
either case, this implicitly limits the source type to being volumetric, as
even in the point source case the source will be "smeared" throughout the
source region that contains the point source coordinate. A source domain is
specified via the ``domains`` constraint placed on the
:class:`openmc.IndependentSource` Python class.
- The source must be isotropic (default for a source)
- The source must use a discrete (i.e., multigroup) energy distribution. The

View file

@ -131,6 +131,12 @@ public:
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
source_region_map_;
// 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_;
// If transport corrected MGXS data is being used, there may be negative
// in-group scattering cross sections that can result in instability in MOC
// and random ray if used naively. This flag enables a stabilization
@ -141,7 +147,7 @@ protected:
//----------------------------------------------------------------------------
// Methods
void apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, int64_t sr);
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);

View file

@ -946,17 +946,16 @@ void FlatSourceDomain::output_to_vtk() const
}
void FlatSourceDomain::apply_external_source_to_source_region(
Discrete* discrete, double strength_factor, int64_t sr)
Discrete* discrete, double strength_factor, SourceRegionHandle& srh)
{
source_regions_.external_source_present(sr) = 1;
srh.external_source_present() = 1;
const auto& discrete_energies = discrete->x();
const auto& discrete_probs = discrete->prob();
for (int i = 0; i < discrete_energies.size(); i++) {
int g = data::mg.get_group_index(discrete_energies[i]);
source_regions_.external_source(sr, g) +=
discrete_probs[i] * strength_factor;
srh.external_source(g) += discrete_probs[i] * strength_factor;
}
}
@ -980,8 +979,9 @@ 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;
apply_external_source_to_source_region(
discrete, strength_factor, source_region);
SourceRegionHandle srh =
source_regions_.get_source_region_handle(source_region);
apply_external_source_to_source_region(discrete, strength_factor, srh);
}
}
}
@ -1023,34 +1023,88 @@ void FlatSourceDomain::convert_external_sources()
{
// Loop over external sources
for (int es = 0; es < model::external_sources.size(); es++) {
// Extract source information
Source* s = model::external_sources[es].get();
IndependentSource* is = dynamic_cast<IndependentSource*>(s);
Discrete* energy = dynamic_cast<Discrete*>(is->energy());
const std::unordered_set<int32_t>& domain_ids = is->domain_ids();
double strength_factor = is->strength();
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);
// If there is no domain constraint specified, then this must be a point
// source. In this case, we need to find the source region that contains the
// point source and apply or relate it to the external source.
if (is->domain_ids().size() == 0) {
// Extract the point source coordinate and find the base source region at
// that point
auto sp = dynamic_cast<SpatialPoint*>(is->space());
GeometryState gs;
gs.r() = sp->r();
gs.r_last() = sp->r();
gs.u() = {1.0, 0.0, 0.0};
bool found = exhaustive_find_cell(gs);
if (!found) {
fatal_error(fmt::format("Could not find cell containing external "
"point source at {}",
sp->r()));
}
int i_cell = gs.lowest_coord().cell;
int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
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);
}
} 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);
}
} 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_) {
} else {
// If not a point source, then use the volumetric domain constraints to
// determine which source regions to apply the external source to.
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);
}
}
} 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);
}
} 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);
}
}
}
}
} // End loop over external sources
@ -1399,6 +1453,25 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
sr_key, {base_source_regions_.get_source_region_handle(sr), sr});
discovered_source_regions_.unlock(sr_key);
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) {
for (int g = 0; g < negroups_; g++) {
double sigma_t = sigma_t_[material * negroups_ + g];
handle.external_source(g) /= sigma_t;
}
}
}
return handle;
}

View file

@ -281,10 +281,21 @@ void validate_random_ray_inputs()
"allowed in random ray mode.");
}
// Validate that a domain ID was specified
if (is->domain_ids().size() == 0) {
fatal_error("Fixed sources must be specified by domain "
"id (cell, material, or universe) in random ray mode.");
// Validate that a domain ID was specified OR that it is a point source
auto sp = dynamic_cast<SpatialPoint*>(is->space());
if (is->domain_ids().size() == 0 && !sp) {
fatal_error("Fixed 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 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.");
}
}
// Check that a discrete energy distribution was used
@ -393,12 +404,12 @@ RandomRaySimulation::RandomRaySimulation()
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_->count_external_source_regions();
}
domain_->apply_meshes();
}
void RandomRaySimulation::prepare_fixed_sources_adjoint(
@ -517,6 +528,8 @@ void RandomRaySimulation::simulate()
finalize_generation();
finalize_batch();
} // End random ray power iteration loop
domain_->count_external_source_regions();
}
void RandomRaySimulation::output_simulation_results() const

View file

@ -0,0 +1,253 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="source">
<density units="macro" value="1.0"/>
<macroscopic name="source"/>
</material>
<material id="2" name="void">
<density units="macro" value="1.0"/>
<macroscopic name="void"/>
</material>
<material id="3" name="absorber">
<density units="macro" value="1.0"/>
<macroscopic name="absorber"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" name="infinite source region" universe="1"/>
<cell id="2" material="2" name="infinite void region" universe="2"/>
<cell id="3" material="3" name="infinite absorber region" universe="3"/>
<cell fill="4" id="4" universe="5"/>
<cell fill="5" id="5" name="full domain" 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 boundary="reflective" coeffs="0.0" id="1" type="x-plane"/>
<surface boundary="vacuum" coeffs="30.0" id="2" type="x-plane"/>
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane"/>
<surface boundary="vacuum" coeffs="30.0" id="4" type="y-plane"/>
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane"/>
<surface boundary="vacuum" coeffs="30.0" id="6" type="z-plane"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>90</particles>
<batches>30</batches>
<inactive>15</inactive>
<source particle="neutron" strength="3.14" type="independent">
<space type="point">
<parameters>2.5 2.5 2.5</parameters>
</space>
<energy type="discrete">
<parameters>100.0 1.0</parameters>
</energy>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>500.0</distance_active>
<distance_inactive>100.0</distance_inactive>
<source particle="neutron" strength="1.0" type="independent">
<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>
<source_region_meshes>
<mesh id="1">
<domain id="6" type="universe"/>
</mesh>
</source_region_meshes>
</random_ray>
<mesh id="1">
<dimension>30 30 30</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
</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>

View file

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

View file

@ -0,0 +1,44 @@
import os
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_point_source_locator():
model = random_ray_three_region_cube()
# Overlay subdivided SR mesh to reduce resolution from 2.5cm -> 1cm
width = 30.0
mesh = openmc.RegularMesh()
mesh.dimension = (30, 30, 30)
mesh.lower_left = (0.0, 0.0, 0.0)
mesh.upper_right = (width, width, width)
model.settings.random_ray['source_region_meshes'] = [
(mesh, [model.geometry.root_universe]),
]
# Define a point source
strengths = [1.0]
midpoints = [100.0]
energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
spatial_distribution = openmc.stats.Point([2.5, 2.5, 2.5])
source = openmc.IndependentSource(
energy=energy_distribution, space=spatial_distribution, strength=3.14)
model.settings.source = [source]
# Settings
model.settings.inactive = 15
model.settings.batches = 30
harness = MGXSTestHarness('statepoint.30.h5', model)
harness.main()