mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Random Ray Base Source Region Refactor (#3576)
This commit is contained in:
John Tramm
GitHub
parent
feefcc6713
commit
3ac64d9a01
17 changed files with 700 additions and 446 deletions
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
//----------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -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};
|
||||
|
||||
|
|
|
|||
|
|
@ -308,7 +308,6 @@ public:
|
|||
//----------------------------------------------------------------------------
|
||||
// Constructors
|
||||
SourceRegion(int negroups, bool is_linear);
|
||||
SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
|
||||
SourceRegion() = default;
|
||||
|
||||
//----------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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(
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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(),
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
244
tests/regression_tests/random_ray_low_density/inputs_true.dat
Normal file
244
tests/regression_tests/random_ray_low_density/inputs_true.dat
Normal 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>
|
||||
|
|
@ -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
60
tests/regression_tests/random_ray_low_density/test.py
Normal 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()
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue