crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 06:25:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Random Ray Source Region Mesh Subdivision (Cell-Under-Voxel Geometry) (#3333)

Browse source 
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
John Tramm 2025-03-06 20:48:31 -06:00 • committed by GitHub
parent e8c9134ff6
commit 9b5678b5f0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
35 changed files with 16638 additions and 460 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
docs/source/_images/2x2_sr_mesh.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 107 KiB

18
docs/source/io_formats/settings.rst
View file

@ -470,6 +470,24 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: prng
:source_region_meshes:
Relates meshes to spatial domains for subdividing source regions with each domain.
:mesh:
Contains an ``id`` attribute and one or more ``<domain>`` sub-elements.
:id:
The unique identifier for the mesh.
:domain:
Each domain element has an ``id`` attribute and a ``type`` attribute.
:id:
The unique identifier for the domain.
:type:
The type of the domain. Can be ``material``, ``cell``, or ``universe``.
----------------------------------
``<resonance_scattering>`` Element
----------------------------------

100
docs/source/usersguide/random_ray.rst
View file

@ -319,20 +319,24 @@ Default behavior using OpenMC's native PRNG can be manually specified as::
.. _subdivision_fsr:
----------------------------------
Subdivision of Flat Source Regions
----------------------------------
-----------------------------
Subdivision of Source Regions
-----------------------------
While the scattering and fission sources in Monte Carlo
are treated continuously, they are assumed to be invariant (flat) within a
MOC or random ray flat source region (FSR). This introduces bias into the
simulation, which can be remedied by reducing the physical size of the FSR
to dimensions below that of typical mean free paths of particles.
While the scattering and fission sources in Monte Carlo are treated
continuously, they are assumed to have a shape (flat or linear) within a MOC or
random ray source region (SR). This introduces bias into the simulation that can
be remedied by reducing the physical size of the SR to be smaller than the
typical mean free paths of particles. While use of linear sources in OpenMC
greatly reduces the error stemming from this approximation, subdivision is still
typically required.
In OpenMC, this subdivision currently must be done manually. The level of
In OpenMC, this subdivision can be done either manually by the user (by defining
additional surfaces and cells in the geometry) or automatically by assigning a
mesh to one or more cells, universes, or material types. The level of
subdivision needed will be dependent on the fidelity the user requires. For
typical light water reactor analysis, consider the following example subdivision
of a two-dimensional 2x2 reflective pincell lattice:
typical light water reactor analysis, consider the following example of manual
subdivision of a two-dimensional 2x2 reflective pincell lattice:
.. figure:: ../_images/2x2_materials.jpeg
:class: with-border
@ -344,9 +348,79 @@ of a two-dimensional 2x2 reflective pincell lattice:
:class: with-border
:width: 400
FSR decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
Manual decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
In the future, automated subdivision of FSRs via mesh overlay may be supported.
Geometry cells can also be subdivided into small source regions by assigning a
mesh to a list of domains, with each domain being of type
:class:`openmc.Material`, :class:`openmc.Cell`, or :class:`openmc.Universe`. The
idea of defining a source region as a combination of a base geometry cell and a
mesh element is known as "cell-under-voxel" style geometry, although in OpenMC
the mesh can be any kind and is not restricted to 3D regular voxels. An example
of overlaying a simple 2D mesh over a geometry is given as::
sr_mesh = openmc.RegularMesh()
sr_mesh.dimension = (n, n)
sr_mesh.lower_left = (0.0, 0.0)
sr_mesh.upper_right = (x, y)
domain = geometry.root_universe
settings.random_ray['source_region_meshes'] = [(sr_mesh, [domain])]
In the above example, we apply a single :math:`n \times n` uniform mesh over the
entire domain by assigning it to the root universe of the geometry.
Alternatively, we might want to apply a finer or coarser mesh to different
regions of a 3D problem, for instance, as::
fuel = openmc.Material(name='UO2 fuel')
...
water = openmc.Material(name='hot borated water')
...
clad = openmc.Material(name='Zr cladding')
...
coarse_mesh = openmc.RegularMesh()
coarse_mesh.dimension = (n, n, n)
coarse_mesh.lower_left = (0.0, 0.0, 0.0)
coarse_mesh.upper_right = (x, y, z)
fine_mesh = openmc.RegularMesh()
fine_mesh.dimension = (2*n, 2*n, 2*n)
fine_mesh.lower_left = (0.0, 0.0, 0.0)
fine_mesh.upper_right = (x, y, z)
settings.random_ray['source_region_meshes'] = [(fine_mesh, [fuel, clad]), (coarse_mesh, [water])]
Note that we don't need to adjust the outer bounds of the mesh to tightly wrap
the domain we assign the mesh to. Rather, OpenMC will dynamically generate
source regions based on the mesh bins rays actually visit, such that no
additional memory is wasted even if a domain only intersects a few mesh bins.
Going back to our 2x2 lattice example, if using a mesh-based subdivision, this
might look as below:
.. figure:: ../_images/2x2_sr_mesh.png
:class: with-border
:width: 400
20x20 overlaid "cell-under-voxel" mesh decomposition for an asymmetrical 2x2 lattice (1.26 cm pitch)
Note that mesh-bashed subdivision is much easier for a user to implement but
does have a few downsides compared to manual subdivision. Manual subdivision can
be done with the specifics of the geometry in mind. As in the pincell example,
it is more efficient to subdivide the fuel region into azimuthal sectors and
radial rings as opposed to a Cartesian mesh. This is more efficient because the
regions are a more uniform size and follow the material boundaries closer,
resulting in the need for fewer source regions. Fewer source regions tends to
equate to a faster computational speed and/or the need for fewer rays per batch
to achieve good statistics. Additionally, applying a mesh often tends to create
a few very small source regions, as shown in the above picture where corners of
the mesh happen to intersect close to the actual fuel-moderator interface. These
small regions are rarely visited by rays, which can result in inaccurate
estimates of the source within those small regions and, thereby, numerical
instability. However, OpenMC utilizes several techniques to detect these small
source regions and mitigate instabilities that are associated with them. In
conclusion, mesh overlay is a great way to subdivide any geometry into smaller
source regions. It can be used while retaining stability, though typically at
the cost of generating more source regions relative to an optimal manual
subdivision.
.. _usersguide_flux_norm:

12
docs/source/usersguide/variance_reduction.rst
View file

@ -130,9 +130,15 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
.. warning::
If using FW-CADIS weight window generation, ensure that the selected weight
window mesh does not subdivide any cells in the problem. In the future, this
restriction is intended to be relaxed, but for now subdivision of cells by a
mesh tally will result in undefined behavior.
window mesh does not subdivide any source regions in the problem. This can
be ensured by assigning the weight window tally mesh to the root universe so
as to create source region boundaries that conform to the mesh, as in the
example below.
::
root = model.geometry.root_universe
settings.random_ray['source_region_meshes'] = [(ww_mesh, [root])]
6. When running your multigroup random ray input deck, OpenMC will automatically
run a forward solve followed by an adjoint solve, with a

4
include/openmc/constants.h
View file

@ -59,6 +59,10 @@ constexpr double RADIAL_MESH_TOL {1e-10};
// Maximum number of random samples per history
constexpr int MAX_SAMPLE {100000};
// Avg. number of hits per batch to be defined as a "small"
// source region in the random ray solver
constexpr double MIN_HITS_PER_BATCH {1.5};
// ============================================================================
// MATH AND PHYSICAL CONSTANTS

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

@ -4,8 +4,10 @@
#include "openmc/constants.h"
#include "openmc/openmp_interface.h"
#include "openmc/position.h"
#include "openmc/random_ray/parallel_map.h"
#include "openmc/random_ray/source_region.h"
#include "openmc/source.h"
#include <unordered_map>
#include <unordered_set>
namespace openmc {
@ -37,7 +39,6 @@ public:
void random_ray_tally();
virtual void accumulate_iteration_flux();
void output_to_vtk() const;
void all_reduce_replicated_source_regions();
void convert_external_sources();
void count_external_source_regions();
void set_adjoint_sources(const vector<double>& forward_flux);
@ -47,12 +48,34 @@ public:
void flatten_xs();
void transpose_scattering_matrix();
void serialize_final_fluxes(vector<double>& flux);
void apply_meshes();
void apply_mesh_to_cell_instances(int32_t i_cell, int32_t mesh_idx,
int target_material_id, const vector<int32_t>& instances,
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);
void finalize_discovered_source_regions();
int64_t n_source_regions() const
{
return source_regions_.n_source_regions();
}
int64_t n_source_elements() const
{
return source_regions_.n_source_regions() * negroups_;
}
//----------------------------------------------------------------------------
// Static Data members
static bool volume_normalized_flux_tallies_;
static bool adjoint_; // If the user wants outputs based on the adjoint flux
// Static variables to store source region meshes and domains
static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
mesh_domain_map_;
//----------------------------------------------------------------------------
// Static data members
static RandomRayVolumeEstimator volume_estimator_;
@ -61,7 +84,6 @@ public:
// Public Data members
bool mapped_all_tallies_ {false}; // If all source regions have been visited
int64_t n_source_regions_ {0}; // Total number of source regions in the model
int64_t n_external_source_regions_ {0}; // Total number of source regions with
// non-zero external source terms
@ -84,6 +106,27 @@ 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.
// It is keyed with a SourceRegionKey, which combines the base source
// region index and the mesh bin.
ParallelMap<SourceRegionKey, SourceRegion, SourceRegionKey::HashFunctor>
discovered_source_regions_;
// Map that relates a SourceRegionKey to the index at which the source
// region can be found in the "source_regions_" container.
std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
source_region_map_;
protected:
//----------------------------------------------------------------------------
// Methods
@ -100,9 +143,7 @@ protected:
//----------------------------------------------------------------------------
// Private data members
int negroups_; // Number of energy groups in simulation
int64_t n_source_elements_ {0}; // Total number of source regions in the model
// times the number of energy groups
int negroups_; // Number of energy groups in simulation
double
simulation_volume_; // Total physical volume of the simulation domain, as

193
include/openmc/random_ray/parallel_map.h Normal file
View file

@ -0,0 +1,193 @@
#ifndef OPENMC_RANDOM_RAY_PARALLEL_HASH_MAP_H
#define OPENMC_RANDOM_RAY_PARALLEL_HASH_MAP_H
#include "openmc/openmp_interface.h"
#include <memory>
#include <unordered_map>
namespace openmc {
/*
* The ParallelMap class allows for threadsafe access to a map-like data
* structure. It is implemented as a hash table with a fixed number of buckets,
* each of which contains a mutex lock and an unordered_map. The class provides
* a subset of the functionality of std::unordered_map. Users must first lock
* the object (using the key) before accessing or modifying the map. The object
* is locked by bucket, allowing for multiple threads to manipulate different
* keys simultaneously, though sometimes threads will need to wait if keys
* happen to be in the same bucket. The ParallelMap will generate pointers to
* hold values, rather than direct storage of values, so as to allow for
* pointers to values to remain valid even after the lock has been released
* (though locking of those values is then left to the user). Iterators to the
* class are provided but are not threadsafe, and are meant to be used only in a
* serial context (e.g., to dump the contents of the map to another data
* structure).
*/
template<typename KeyType, typename ValueType, typename HashFunctor>
class ParallelMap {
//----------------------------------------------------------------------------
// Helper structs and classes
struct Bucket {
OpenMPMutex lock_;
std::unordered_map<KeyType, std::unique_ptr<ValueType>, HashFunctor> map_;
};
// The iterator yields a pair: (const KeyType&, ValueType&)
class iterator {
public:
using iterator_category = std::forward_iterator_tag;
using value_type = std::pair<const KeyType&, ValueType&>;
using difference_type = std::ptrdiff_t;
using pointer = void; // Not providing pointer semantics.
using reference = value_type;
iterator(std::vector<Bucket>* buckets, std::size_t bucket_index,
typename std::unordered_map<KeyType, std::unique_ptr<ValueType>,
HashFunctor>::iterator inner_it)
: buckets_(buckets), bucket_index_(bucket_index), inner_it_(inner_it)
{
// Advance to the first valid element if necessary.
advance_to_valid();
}
// Dereference returns a pair of (key, value).
reference operator*() const
{
return {inner_it_->first, *inner_it_->second};
}
iterator& operator++()
{
++inner_it_;
advance_to_valid();
return *this;
}
iterator operator++(int)
{
iterator tmp = *this;
++(*this);
return tmp;
}
bool operator==(const iterator& other) const
{
// Two iterators are equal if they refer to the same bucket vector and are
// both at end, or if they have the same bucket index and inner iterator.
return buckets_ == other.buckets_ &&
bucket_index_ == other.bucket_index_ &&
(bucket_index_ == buckets_->size() ||
inner_it_ == other.inner_it_);
}
bool operator!=(const iterator& other) const { return !(*this == other); }
private:
// Helper function: if we are at the end of the current bucket, advance to
// the next non-empty bucket.
void advance_to_valid()
{
while (bucket_index_ < buckets_->size() &&
inner_it_ == (*buckets_)[bucket_index_].map_.end()) {
++bucket_index_;
if (bucket_index_ < buckets_->size())
inner_it_ = (*buckets_)[bucket_index_].map_.begin();
}
}
std::vector<Bucket>* buckets_;
std::size_t bucket_index_;
typename std::unordered_map<KeyType, std::unique_ptr<ValueType>,
HashFunctor>::iterator inner_it_;
};
public:
//----------------------------------------------------------------------------
// Constructor
ParallelMap(int n_buckets = 1000) : buckets_(n_buckets) {}
//----------------------------------------------------------------------------
// Public Methods
void lock(const KeyType& key)
{
Bucket& bucket = get_bucket(key);
bucket.lock_.lock();
}
void unlock(const KeyType& key)
{
Bucket& bucket = get_bucket(key);
bucket.lock_.unlock();
}
void clear()
{
for (auto& bucket : buckets_) {
bucket.map_.clear();
}
}
bool contains(const KeyType& key)
{
Bucket& bucket = get_bucket(key);
// C++20
// return bucket.map_.contains(key);
return bucket.map_.find(key) != bucket.map_.end();
}
ValueType& operator[](const KeyType& key)
{
Bucket& bucket = get_bucket(key);
return *bucket.map_[key].get();
}
ValueType* emplace(KeyType key, const ValueType& value)
{
Bucket& bucket = get_bucket(key);
// Attempt to emplace the new element into the unordered_map within the
auto result = bucket.map_.emplace(key, std::make_unique<ValueType>(value));
auto it = result.first;
return it->second.get();
}
// Return iterator to first element.
iterator begin()
{
std::size_t bucket_index = 0;
auto inner_it = buckets_.empty()
? typename std::unordered_map<KeyType,
std::unique_ptr<ValueType>, HashFunctor>::iterator()
: buckets_[0].map_.begin();
return iterator(&buckets_, bucket_index, inner_it);
}
// Return iterator to one-past-last element.
iterator end()
{
// End is signaled by bucket_index_ equal to buckets_.size()
return iterator(&buckets_, buckets_.size(),
typename std::unordered_map<KeyType, std::unique_ptr<ValueType>,
HashFunctor>::iterator());
}
private:
//----------------------------------------------------------------------------
// Private Methods
Bucket& get_bucket(const KeyType& key)
{
return buckets_[hash(key) % buckets_.size()];
}
//----------------------------------------------------------------------------
// Private Data Fields
HashFunctor hash;
vector<Bucket> buckets_;
};
} // namespace openmc
#endif // OPENMC_RANDOM_RAY_PARALLEL_HASH_MAP_H

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

@ -25,11 +25,17 @@ public:
//----------------------------------------------------------------------------
// Methods
void event_advance_ray();
void attenuate_flux(double distance, bool is_active);
void attenuate_flux_flat_source(double distance, bool is_active);
void attenuate_flux_flat_source_void(double distance, bool is_active);
void attenuate_flux_linear_source(double distance, bool is_active);
void attenuate_flux_linear_source_void(double distance, bool is_active);
void attenuate_flux(double distance, bool is_active, double offset = 0.0);
void attenuate_flux_inner(
double distance, bool is_active, int64_t sr, int mesh_bin, Position r);
void attenuate_flux_flat_source(
SourceRegionHandle& srh, double distance, bool is_active, Position r);
void attenuate_flux_flat_source_void(
SourceRegionHandle& srh, double distance, bool is_active, Position r);
void attenuate_flux_linear_source(
SourceRegionHandle& srh, double distance, bool is_active, Position r);
void attenuate_flux_linear_source_void(
SourceRegionHandle& srh, double distance, bool is_active, Position r);
void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
uint64_t transport_history_based_single_ray();
@ -42,6 +48,7 @@ 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
//----------------------------------------------------------------------------
@ -55,6 +62,8 @@ private:
// Private data members
vector<float> delta_psi_;
vector<MomentArray> delta_moments_;
vector<int> mesh_bins_;
vector<double> mesh_fractional_lengths_;
int negroups_;
FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source

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

@ -20,10 +20,13 @@ public:
//----------------------------------------------------------------------------
// Methods
void compute_segment_correction_factors();
void prepare_fixed_sources();
void prepare_fixed_sources_adjoint(vector<double>& forward_flux);
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 simulate();
void reduce_simulation_statistics();
void output_simulation_results() const;
void instability_check(
int64_t n_hits, double k_eff, double& avg_miss_rate) const;
@ -63,6 +66,7 @@ private:
void openmc_run_random_ray();
void validate_random_ray_inputs();
void openmc_reset_random_ray();
} // namespace openmc

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

@ -44,7 +44,7 @@ inline void hash_combine(size_t& seed, const size_t v)
}
//----------------------------------------------------------------------------
// Helper Structs
// Helper Structs and Classes
// A mapping object that is used to map between a specific random ray
// source region and an OpenMC native tally bin that it should score to
@ -81,11 +81,234 @@ struct TallyTask {
};
};
// The SourceRegionKey combines a base source region (i.e., a material
// filled cell instance) with a mesh bin. This key is used as a handle
// for dynamically discovered source regions when subdividing source
// regions with meshes.
class SourceRegionKey {
public:
int64_t base_source_region_id;
int64_t mesh_bin;
SourceRegionKey() = default;
SourceRegionKey(int64_t source_region, int64_t bin)
: base_source_region_id(source_region), mesh_bin(bin)
{}
// Equality operator required by the unordered_map
bool operator==(const SourceRegionKey& other) const
{
return base_source_region_id == other.base_source_region_id &&
mesh_bin == other.mesh_bin;
}
// Less than operator required by std::sort
bool operator<(const SourceRegionKey& other) const
{
if (base_source_region_id < other.base_source_region_id) {
return true;
} else if (base_source_region_id > other.base_source_region_id) {
return false;
} else {
return mesh_bin < other.mesh_bin;
}
}
// Hashing functor required by the unordered_map
struct HashFunctor {
size_t operator()(const SourceRegionKey& key) const
{
size_t seed = 0;
hash_combine(seed, key.base_source_region_id);
hash_combine(seed, key.mesh_bin);
return seed;
}
};
};
// Forward declaration of SourceRegion
class SourceRegion;
class SourceRegionHandle {
public:
//----------------------------------------------------------------------------
// Constructors
SourceRegionHandle(SourceRegion& sr);
SourceRegionHandle() = default;
// All fields are commented/described in the SourceRegion class definition
// below
//----------------------------------------------------------------------------
// Public Data members
int negroups_;
bool is_numerical_fp_artifact_ {false};
bool is_linear_ {false};
// Scalar fields
int* material_;
int* is_small_;
int* n_hits_;
int* birthday_;
OpenMPMutex* lock_;
double* volume_;
double* volume_t_;
double* volume_sq_;
double* volume_sq_t_;
double* volume_naive_;
int* position_recorded_;
int* external_source_present_;
Position* position_;
Position* centroid_;
Position* centroid_iteration_;
Position* centroid_t_;
MomentMatrix* mom_matrix_;
MomentMatrix* mom_matrix_t_;
// A set of volume tally tasks. This more complicated data structure is
// convenient for ensuring that volumes are only tallied once per source
// region, regardless of how many energy groups are used for tallying.
std::unordered_set<TallyTask, TallyTask::HashFunctor>* volume_task_;
// Mesh that subdivides this source region
int* mesh_;
int64_t* parent_sr_;
// Energy group-wise 1D arrays
double* scalar_flux_old_;
double* scalar_flux_new_;
float* source_;
float* external_source_;
double* scalar_flux_final_;
MomentArray* source_gradients_;
MomentArray* flux_moments_old_;
MomentArray* flux_moments_new_;
MomentArray* flux_moments_t_;
// 2D array representing values for all energy groups x tally
// tasks. Each group may have a different number of tally tasks
// associated with it, necessitating the use of a jagged array.
vector<TallyTask>* tally_task_;
//----------------------------------------------------------------------------
// Public Accessors
int& material() { return *material_; }
const int material() const { return *material_; }
int& is_small() { return *is_small_; }
const int is_small() const { return *is_small_; }
int& n_hits() { return *n_hits_; }
const int n_hits() const { return *n_hits_; }
void lock() { lock_->lock(); }
void unlock() { lock_->unlock(); }
double& volume() { return *volume_; }
const double volume() const { return *volume_; }
double& volume_t() { return *volume_t_; }
const double volume_t() const { return *volume_t_; }
double& volume_sq() { return *volume_sq_; }
const double volume_sq() const { return *volume_sq_; }
double& volume_sq_t() { return *volume_sq_t_; }
const double volume_sq_t() const { return *volume_sq_t_; }
double& volume_naive() { return *volume_naive_; }
const double volume_naive() const { return *volume_naive_; }
int& position_recorded() { return *position_recorded_; }
const int position_recorded() const { return *position_recorded_; }
int& external_source_present() { return *external_source_present_; }
const int external_source_present() const
{
return *external_source_present_;
}
Position& position() { return *position_; }
const Position position() const { return *position_; }
Position& centroid() { return *centroid_; }
const Position centroid() const { return *centroid_; }
Position& centroid_iteration() { return *centroid_iteration_; }
const Position centroid_iteration() const { return *centroid_iteration_; }
Position& centroid_t() { return *centroid_t_; }
const Position centroid_t() const { return *centroid_t_; }
MomentMatrix& mom_matrix() { return *mom_matrix_; }
const MomentMatrix mom_matrix() const { return *mom_matrix_; }
MomentMatrix& mom_matrix_t() { return *mom_matrix_t_; }
const MomentMatrix mom_matrix_t() const { return *mom_matrix_t_; }
std::unordered_set<TallyTask, TallyTask::HashFunctor>& volume_task()
{
return *volume_task_;
}
const std::unordered_set<TallyTask, TallyTask::HashFunctor>& volume_task()
const
{
return *volume_task_;
}
int& mesh() { return *mesh_; }
const int mesh() const { return *mesh_; }
int64_t& parent_sr() { return *parent_sr_; }
const int64_t parent_sr() const { return *parent_sr_; }
double& scalar_flux_old(int g) { return scalar_flux_old_[g]; }
const double scalar_flux_old(int g) const { return scalar_flux_old_[g]; }
double& scalar_flux_new(int g) { return scalar_flux_new_[g]; }
const double scalar_flux_new(int g) const { return scalar_flux_new_[g]; }
double& scalar_flux_final(int g) { return scalar_flux_final_[g]; }
const double scalar_flux_final(int g) const { return scalar_flux_final_[g]; }
float& source(int g) { return source_[g]; }
const float source(int g) const { return source_[g]; }
float& external_source(int g) { return external_source_[g]; }
const float external_source(int g) const { return external_source_[g]; }
MomentArray& source_gradients(int g) { return source_gradients_[g]; }
const MomentArray source_gradients(int g) const
{
return source_gradients_[g];
}
MomentArray& flux_moments_old(int g) { return flux_moments_old_[g]; }
const MomentArray flux_moments_old(int g) const
{
return flux_moments_old_[g];
}
MomentArray& flux_moments_new(int g) { return flux_moments_new_[g]; }
const MomentArray flux_moments_new(int g) const
{
return flux_moments_new_[g];
}
MomentArray& flux_moments_t(int g) { return flux_moments_t_[g]; }
const MomentArray flux_moments_t(int g) const { return flux_moments_t_[g]; }
vector<TallyTask>& tally_task(int g) { return tally_task_[g]; }
const vector<TallyTask>& tally_task(int g) const { return tally_task_[g]; }
}; // class SourceRegionHandle
class SourceRegion {
public:
//----------------------------------------------------------------------------
// Constructors
SourceRegion(int negroups, bool is_linear);
SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr);
SourceRegion() = default;
//----------------------------------------------------------------------------
@ -104,7 +327,13 @@ public:
double volume_naive_ {0.0}; //!< Volume as integrated from this iteration only
int position_recorded_ {0}; //!< Has the position been recorded yet?
int external_source_present_ {
0}; //!< Is an external source present in this region?
0}; //!< Is an external source present in this region?
int is_small_ {0}; //!< Is it "small", receiving < 1.5 hits per iteration?
int n_hits_ {0}; //!< Number of total hits (ray crossings)
// Mesh that subdivides this source region
int mesh_ {C_NONE}; //!< Index in openmc::model::meshes array that subdivides
//!< this source region
int64_t parent_sr_ {C_NONE}; //!< Index of a parent source region
Position position_ {
0.0, 0.0, 0.0}; //!< A position somewhere inside the region
Position centroid_ {0.0, 0.0, 0.0}; //!< The centroid
@ -150,7 +379,6 @@ public:
// tasks. Each group may have a different number of tally tasks
// associated with it, necessitating the use of a jagged array.
vector<vector<TallyTask>> tally_task_;
}; // class SourceRegion
class SourceRegionContainer {
@ -165,28 +393,34 @@ public:
//----------------------------------------------------------------------------
// Public Accessors
int& material(int64_t sr) { return material_[sr]; }
const int& material(int64_t sr) const { return material_[sr]; }
const int material(int64_t sr) const { return material_[sr]; }
int& is_small(int64_t sr) { return is_small_[sr]; }
const int is_small(int64_t sr) const { return is_small_[sr]; }
int& n_hits(int64_t sr) { return n_hits_[sr]; }
const int n_hits(int64_t sr) const { return n_hits_[sr]; }
OpenMPMutex& lock(int64_t sr) { return lock_[sr]; }
const OpenMPMutex& lock(int64_t sr) const { return lock_[sr]; }
double& volume(int64_t sr) { return volume_[sr]; }
const double& volume(int64_t sr) const { return volume_[sr]; }
const double volume(int64_t sr) const { return volume_[sr]; }
double& volume_t(int64_t sr) { return volume_t_[sr]; }
const double& volume_t(int64_t sr) const { return volume_t_[sr]; }
const double volume_t(int64_t sr) const { return volume_t_[sr]; }
double& volume_sq(int64_t sr) { return volume_sq_[sr]; }
const double& volume_sq(int64_t sr) const { return volume_sq_[sr]; }
const double volume_sq(int64_t sr) const { return volume_sq_[sr]; }
double& volume_sq_t(int64_t sr) { return volume_sq_t_[sr]; }
const double& volume_sq_t(int64_t sr) const { return volume_sq_t_[sr]; }
const double volume_sq_t(int64_t sr) const { return volume_sq_t_[sr]; }
double& volume_naive(int64_t sr) { return volume_naive_[sr]; }
const double& volume_naive(int64_t sr) const { return volume_naive_[sr]; }
const double volume_naive(int64_t sr) const { return volume_naive_[sr]; }
int& position_recorded(int64_t sr) { return position_recorded_[sr]; }
const int& position_recorded(int64_t sr) const
const int position_recorded(int64_t sr) const
{
return position_recorded_[sr];
}
@ -195,31 +429,31 @@ public:
{
return external_source_present_[sr];
}
const int& external_source_present(int64_t sr) const
const int external_source_present(int64_t sr) const
{
return external_source_present_[sr];
}
Position& position(int64_t sr) { return position_[sr]; }
const Position& position(int64_t sr) const { return position_[sr]; }
const Position position(int64_t sr) const { return position_[sr]; }
Position& centroid(int64_t sr) { return centroid_[sr]; }
const Position& centroid(int64_t sr) const { return centroid_[sr]; }
const Position centroid(int64_t sr) const { return centroid_[sr]; }
Position& centroid_iteration(int64_t sr) { return centroid_iteration_[sr]; }
const Position& centroid_iteration(int64_t sr) const
const Position centroid_iteration(int64_t sr) const
{
return centroid_iteration_[sr];
}
Position& centroid_t(int64_t sr) { return centroid_t_[sr]; }
const Position& centroid_t(int64_t sr) const { return centroid_t_[sr]; }
const Position centroid_t(int64_t sr) const { return centroid_t_[sr]; }
MomentMatrix& mom_matrix(int64_t sr) { return mom_matrix_[sr]; }
const MomentMatrix& mom_matrix(int64_t sr) const { return mom_matrix_[sr]; }
const MomentMatrix mom_matrix(int64_t sr) const { return mom_matrix_[sr]; }
MomentMatrix& mom_matrix_t(int64_t sr) { return mom_matrix_t_[sr]; }
const MomentMatrix& mom_matrix_t(int64_t sr) const
const MomentMatrix mom_matrix_t(int64_t sr) const
{
return mom_matrix_t_[sr];
}
@ -228,12 +462,12 @@ public:
{
return source_gradients_[index(sr, g)];
}
const MomentArray& source_gradients(int64_t sr, int g) const
const MomentArray source_gradients(int64_t sr, int g) const
{
return source_gradients_[index(sr, g)];
}
MomentArray& source_gradients(int64_t se) { return source_gradients_[se]; }
const MomentArray& source_gradients(int64_t se) const
const MomentArray source_gradients(int64_t se) const
{
return source_gradients_[se];
}
@ -242,12 +476,12 @@ public:
{
return flux_moments_old_[index(sr, g)];
}
const MomentArray& flux_moments_old(int64_t sr, int g) const
const MomentArray flux_moments_old(int64_t sr, int g) const
{
return flux_moments_old_[index(sr, g)];
}
MomentArray& flux_moments_old(int64_t se) { return flux_moments_old_[se]; }
const MomentArray& flux_moments_old(int64_t se) const
const MomentArray flux_moments_old(int64_t se) const
{
return flux_moments_old_[se];
}
@ -256,12 +490,12 @@ public:
{
return flux_moments_new_[index(sr, g)];
}
const MomentArray& flux_moments_new(int64_t sr, int g) const
const MomentArray flux_moments_new(int64_t sr, int g) const
{
return flux_moments_new_[index(sr, g)];
}
MomentArray& flux_moments_new(int64_t se) { return flux_moments_new_[se]; }
const MomentArray& flux_moments_new(int64_t se) const
const MomentArray flux_moments_new(int64_t se) const
{
return flux_moments_new_[se];
}
@ -270,12 +504,12 @@ public:
{
return flux_moments_t_[index(sr, g)];
}
const MomentArray& flux_moments_t(int64_t sr, int g) const
const MomentArray flux_moments_t(int64_t sr, int g) const
{
return flux_moments_t_[index(sr, g)];
}
MomentArray& flux_moments_t(int64_t se) { return flux_moments_t_[se]; }
const MomentArray& flux_moments_t(int64_t se) const
const MomentArray flux_moments_t(int64_t se) const
{
return flux_moments_t_[se];
}
@ -284,12 +518,12 @@ public:
{
return scalar_flux_old_[index(sr, g)];
}
const double& scalar_flux_old(int64_t sr, int g) const
const double scalar_flux_old(int64_t sr, int g) const
{
return scalar_flux_old_[index(sr, g)];
}
double& scalar_flux_old(int64_t se) { return scalar_flux_old_[se]; }
const double& scalar_flux_old(int64_t se) const
const double scalar_flux_old(int64_t se) const
{
return scalar_flux_old_[se];
}
@ -298,12 +532,12 @@ public:
{
return scalar_flux_new_[index(sr, g)];
}
const double& scalar_flux_new(int64_t sr, int g) const
const double scalar_flux_new(int64_t sr, int g) const
{
return scalar_flux_new_[index(sr, g)];
}
double& scalar_flux_new(int64_t se) { return scalar_flux_new_[se]; }
const double& scalar_flux_new(int64_t se) const
const double scalar_flux_new(int64_t se) const
{
return scalar_flux_new_[se];
}
@ -312,34 +546,31 @@ public:
{
return scalar_flux_final_[index(sr, g)];
}
const double& scalar_flux_final(int64_t sr, int g) const
const double scalar_flux_final(int64_t sr, int g) const
{
return scalar_flux_final_[index(sr, g)];
}
double& scalar_flux_final(int64_t se) { return scalar_flux_final_[se]; }
const double& scalar_flux_final(int64_t se) const
const double scalar_flux_final(int64_t se) const
{
return scalar_flux_final_[se];
}
float& source(int64_t sr, int g) { return source_[index(sr, g)]; }
const float& source(int64_t sr, int g) const { return source_[index(sr, g)]; }
const float source(int64_t sr, int g) const { return source_[index(sr, g)]; }
float& source(int64_t se) { return source_[se]; }
const float& source(int64_t se) const { return source_[se]; }
const float source(int64_t se) const { return source_[se]; }
float& external_source(int64_t sr, int g)
{
return external_source_[index(sr, g)];
}
const float& external_source(int64_t sr, int g) const
const float external_source(int64_t sr, int g) const
{
return external_source_[index(sr, g)];
}
float& external_source(int64_t se) { return external_source_[se]; }
const float& external_source(int64_t se) const
{
return external_source_[se];
}
const float external_source(int64_t se) const { return external_source_[se]; }
vector<TallyTask>& tally_task(int64_t sr, int g)
{
@ -365,13 +596,26 @@ public:
return volume_task_[sr];
}
int& mesh(int64_t sr) { return mesh_[sr]; }
const int mesh(int64_t sr) const { return mesh_[sr]; }
int64_t& parent_sr(int64_t sr) { return parent_sr_[sr]; }
const int64_t parent_sr(int64_t sr) const { return parent_sr_[sr]; }
//----------------------------------------------------------------------------
// Public Methods
void push_back(const SourceRegion& sr);
void assign(int n_source_regions, const SourceRegion& source_region);
void flux_swap();
void mpi_sync_ranks(bool reduce_position);
int64_t n_source_regions() const { return n_source_regions_; }
int64_t n_source_elements() const { return n_source_regions_ * negroups_; }
int& negroups() { return negroups_; }
const int negroups() const { return negroups_; }
bool& is_linear() { return is_linear_; }
const bool is_linear() const { return is_linear_; }
SourceRegionHandle get_source_region_handle(int64_t sr);
void adjoint_reset();
private:
//----------------------------------------------------------------------------
@ -382,6 +626,10 @@ private:
// SoA storage for scalar fields (one item per source region)
vector<int> material_;
vector<int> is_small_;
vector<int> n_hits_;
vector<int> mesh_;
vector<int64_t> parent_sr_;
vector<OpenMPMutex> lock_;
vector<double> volume_;
vector<double> volume_t_;

55
openmc/settings.py
View file

@ -7,11 +7,12 @@ from pathlib import Path
import lxml.etree as ET
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.stats.multivariate import MeshSpatial
from ._xml import clean_indentation, get_text, reorder_attributes
from .mesh import _read_meshes, RegularMesh
from .mesh import _read_meshes, RegularMesh, MeshBase
from .source import SourceBase, MeshSource, IndependentSource
from .utility_funcs import input_path
from .volume import VolumeCalculation
@ -180,6 +181,12 @@ class Settings:
:sample_method:
Sampling method for the ray starting location and direction of
travel. Options are `prng` (default) or 'halton`.
:source_region_meshes:
List of tuples where each tuple contains a mesh and a list of
domains. Each domain is an instance of openmc.Material, openmc.Cell,
or openmc.Universe. The mesh will be applied to the listed domains
to subdivide source regions so as to improve accuracy and/or conform
with tally meshes.
.. versionadded:: 0.15.0
resonance_scattering : dict
@ -1156,6 +1163,18 @@ class Settings:
cv.check_type('volume normalized flux tallies', value, bool)
elif key == 'adjoint':
cv.check_type('adjoint', value, bool)
elif key == 'source_region_meshes':
cv.check_type('source region meshes', value, Iterable)
for mesh, domains in value:
cv.check_type('mesh', mesh, MeshBase)
cv.check_type('domains', domains, Iterable)
valid_types = (openmc.Material, openmc.Cell, openmc.Universe)
for domain in domains:
if not isinstance(domain, valid_types):
raise ValueError(
f'Invalid domain type: {type(domain)}. Expected '
'openmc.Material, openmc.Cell, or openmc.Universe.')
cv.check_type('adjoint', value, bool)
elif key == 'sample_method':
cv.check_value('sample method', value,
('prng', 'halton'))
@ -1573,7 +1592,7 @@ class Settings:
root.append(wwg.mesh.to_xml_element())
if mesh_memo is not None:
mesh_memo.add(wwg.mesh)
mesh_memo.add(wwg.mesh.id)
def _create_weight_windows_file_element(self, root):
if self.weight_windows_file is not None:
@ -1604,13 +1623,25 @@ class Settings:
elem = ET.SubElement(root, "max_tracks")
elem.text = str(self._max_tracks)
def _create_random_ray_subelement(self, root):
def _create_random_ray_subelement(self, root, mesh_memo=None):
if self._random_ray:
element = ET.SubElement(root, "random_ray")
for key, value in self._random_ray.items():
if key == 'ray_source' and isinstance(value, SourceBase):
source_element = value.to_xml_element()
element.append(source_element)
elif key == 'source_region_meshes':
subelement = ET.SubElement(element, 'source_region_meshes')
for mesh, domains in value:
mesh_elem = ET.SubElement(subelement, 'mesh')
mesh_elem.set('id', str(mesh.id))
for domain in domains:
domain_elem = ET.SubElement(mesh_elem, 'domain')
domain_elem.set('id', str(domain.id))
domain_elem.set('type', domain.__class__.__name__.lower())
if mesh_memo is not None and mesh.id not in mesh_memo:
root.append(mesh.to_xml_element())
mesh_memo.add(mesh.id)
else:
subelement = ET.SubElement(element, key)
subelement.text = str(value)
@ -2007,6 +2038,22 @@ class Settings:
)
elif child.tag == 'sample_method':
self.random_ray['sample_method'] = child.text
elif child.tag == 'source_region_meshes':
self.random_ray['source_region_meshes'] = []
for mesh_elem in child.findall('mesh'):
mesh = MeshBase.from_xml_element(mesh_elem)
domains = []
for domain_elem in mesh_elem.findall('domain'):
domain_id = int(domain_elem.get('id'))
domain_type = domain_elem.get('type')
if domain_type == 'material':
domain = openmc.Material(domain_id)
elif domain_type == 'cell':
domain = openmc.Cell(domain_id)
elif domain_type == 'universe':
domain = openmc.Universe(domain_id)
domains.append(domain)
self.random_ray['source_region_meshes'].append((mesh, domains))
def _use_decay_photons_from_xml_element(self, root):
text = get_text(root, 'use_decay_photons')
@ -2077,7 +2124,7 @@ class Settings:
self._create_weight_window_checkpoints_subelement(element)
self._create_max_history_splits_subelement(element)
self._create_max_tracks_subelement(element)
self._create_random_ray_subelement(element)
self._create_random_ray_subelement(element, mesh_memo)
self._create_use_decay_photons_subelement(element)
# Clean the indentation in the file to be user-readable

3
src/finalize.cpp
View file

@ -17,6 +17,7 @@
#include "openmc/photon.h"
#include "openmc/plot.h"
#include "openmc/random_lcg.h"
#include "openmc/random_ray/random_ray_simulation.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
@ -174,6 +175,8 @@ int openmc_finalize()
MPI_Type_free(&mpi::source_site);
#endif
openmc_reset_random_ray();
return 0;
}

466
src/random_ray/flat_source_domain.cpp
View file

@ -1,6 +1,7 @@
#include "openmc/random_ray/flat_source_domain.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
@ -14,6 +15,7 @@
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/timer.h"
#include "openmc/weight_windows.h"
#include <cstdio>
@ -28,6 +30,8 @@ RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {
RandomRayVolumeEstimator::HYBRID};
bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
bool FlatSourceDomain::adjoint_ {false};
std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
FlatSourceDomain::mesh_domain_map_;
FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
{
@ -35,20 +39,21 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
// indices, and store the material type The reason for the offsets is that
// some cell types may not have material fills, and therefore do not
// produce FSRs. Thus, we cannot index into the global arrays directly
int base_source_regions = 0;
for (const auto& c : model::cells) {
if (c->type_ != Fill::MATERIAL) {
source_region_offsets_.push_back(-1);
} else {
source_region_offsets_.push_back(n_source_regions_);
n_source_regions_ += c->n_instances_;
n_source_elements_ += c->n_instances_ * negroups_;
source_region_offsets_.push_back(base_source_regions);
base_source_regions += c->n_instances_;
}
}
// Initialize cell-wise arrays
// Initialize source regions.
bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT;
source_regions_ = SourceRegionContainer(negroups_, is_linear);
source_regions_.assign(n_source_regions_, SourceRegion(negroups_, is_linear));
source_regions_.assign(
base_source_regions, SourceRegion(negroups_, is_linear));
// Initialize materials
int64_t source_region_id = 0;
@ -62,7 +67,7 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
}
// Sanity check
if (source_region_id != n_source_regions_) {
if (source_region_id != base_source_regions) {
fatal_error("Unexpected number of source regions");
}
@ -92,12 +97,13 @@ void FlatSourceDomain::batch_reset()
{
// Reset scalar fluxes and iteration volume tallies to zero
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
source_regions_.volume(sr) = 0.0;
source_regions_.volume_sq(sr) = 0.0;
}
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.scalar_flux_new(se) = 0.0;
}
}
@ -105,7 +111,7 @@ void FlatSourceDomain::batch_reset()
void FlatSourceDomain::accumulate_iteration_flux()
{
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.scalar_flux_final(se) +=
source_regions_.scalar_flux_new(se);
}
@ -121,13 +127,13 @@ void FlatSourceDomain::update_neutron_source(double 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++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 0.0;
}
// Add scattering + fission source
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
@ -155,7 +161,7 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
// 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++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
}
}
@ -174,14 +180,14 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes(
// Normalize scalar flux to total distance travelled by all rays this
// iteration
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.scalar_flux_new(se) *= normalization_factor;
}
// Accumulate cell-wise ray length tallies collected this iteration, then
// update the simulation-averaged cell-wise volume estimates
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
source_regions_.volume_t(sr) += source_regions_.volume(sr);
source_regions_.volume_sq_t(sr) += source_regions_.volume_sq(sr);
source_regions_.volume_naive(sr) =
@ -225,9 +231,10 @@ void FlatSourceDomain::set_flux_to_source(int64_t sr, int g)
int64_t FlatSourceDomain::add_source_to_scalar_flux()
{
int64_t n_hits = 0;
double inverse_batch = 1.0 / simulation::current_batch;
#pragma omp parallel for reduction(+ : n_hits)
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
double volume_simulation_avg = source_regions_.volume(sr);
double volume_iteration = source_regions_.volume_naive(sr);
@ -237,6 +244,14 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
n_hits++;
}
// Set the SR to small status if its expected number of hits
// per iteration is less than 1.5
if (source_regions_.n_hits(sr) * inverse_batch < MIN_HITS_PER_BATCH) {
source_regions_.is_small(sr) = 1;
} else {
source_regions_.is_small(sr) = 0;
}
// The volume treatment depends on the volume estimator type
// and whether or not an external source is present in the cell.
double volume;
@ -248,7 +263,8 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
volume = volume_simulation_avg;
break;
case RandomRayVolumeEstimator::HYBRID:
if (source_regions_.external_source_present(sr)) {
if (source_regions_.external_source_present(sr) ||
source_regions_.is_small(sr)) {
volume = volume_iteration;
} else {
volume = volume_simulation_avg;
@ -279,16 +295,12 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux()
// in the cell we will use the previous iteration's flux estimate. This
// injects a small degree of correlation into the simulation, but this
// is going to be trivial when the miss rate is a few percent or less.
if (source_regions_.external_source_present(sr)) {
if (source_regions_.external_source_present(sr) && !adjoint_) {
set_flux_to_old_flux(sr, g);
} else {
set_flux_to_source(sr, g);
}
}
// If the FSR was not hit this iteration, and it has never been hit in
// any iteration (i.e., volume is zero), then we want to set this to 0
// to avoid dividing anything by a zero volume. This happens implicitly
// given that the new scalar flux arrays are set to zero each iteration.
}
}
@ -304,10 +316,10 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
double fission_rate_new = 0;
// Vector for gathering fission source terms for Shannon entropy calculation
vector<float> p(n_source_regions_, 0.0f);
vector<float> p(n_source_regions(), 0.0f);
#pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new)
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
// If simulation averaged volume is zero, don't include this cell
double volume = source_regions_.volume(sr);
@ -350,7 +362,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const
double inverse_sum = 1 / fission_rate_new;
#pragma omp parallel for reduction(+ : H)
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
// Only if FSR has non-negative and non-zero fission source
if (p[sr] > 0.0f) {
// Normalize to total weight of bank sites. p_i for better performance
@ -408,7 +420,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
// Attempt to generate mapping for all source regions
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
// If this source region has not been hit by a ray yet, then
// we aren't going to be able to map it, so skip it.
@ -423,6 +435,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies()
Particle p;
p.r() = source_regions_.position(sr);
p.r_last() = source_regions_.position(sr);
p.u() = {1.0, 0.0, 0.0};
bool found = exhaustive_find_cell(p);
// Loop over energy groups (so as to support energy filters)
@ -517,7 +530,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const
// total external source strength in the simulation.
double simulation_external_source_strength = 0.0;
#pragma omp parallel for reduction(+ : simulation_external_source_strength)
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
double volume = source_regions_.volume(sr) * simulation_volume_;
for (int g = 0; g < negroups_; g++) {
@ -570,7 +583,7 @@ void FlatSourceDomain::random_ray_tally()
// element, we check if there are any scores needed and apply
// them.
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
// The fsr.volume_ is the unitless fractional simulation averaged volume
// (i.e., it is the FSR's fraction of the overall simulation volume). The
// simulation_volume_ is the total 3D physical volume in cm^3 of the
@ -673,30 +686,6 @@ void FlatSourceDomain::random_ray_tally()
openmc::simulation::time_tallies.stop();
}
void FlatSourceDomain::all_reduce_replicated_source_regions()
{
#ifdef OPENMC_MPI
// If we only have 1 MPI rank, no need
// to reduce anything.
if (mpi::n_procs <= 1)
return;
simulation::time_bank_sendrecv.start();
// First, we broadcast the fully mapped tally status variable so that
// all ranks are on the same page
int mapped_all_tallies_i = static_cast<int>(mapped_all_tallies_);
MPI_Bcast(&mapped_all_tallies_i, 1, MPI_INT, 0, mpi::intracomm);
bool reduce_position =
simulation::current_batch > settings::n_inactive && !mapped_all_tallies_i;
source_regions_.mpi_sync_ranks(reduce_position);
simulation::time_bank_sendrecv.stop();
#endif
}
double FlatSourceDomain::evaluate_flux_at_point(
Position r, int64_t sr, int g) const
{
@ -765,8 +754,9 @@ void FlatSourceDomain::output_to_vtk() const
// Relate voxel spatial locations to random ray source regions
vector<int> voxel_indices(Nx * Ny * Nz);
vector<Position> voxel_positions(Nx * Ny * Nz);
#pragma omp parallel for collapse(3)
vector<double> weight_windows(Nx * Ny * Nz);
float min_weight = 1e20;
#pragma omp parallel for collapse(3) reduction(min : min_weight)
for (int z = 0; z < Nz; z++) {
for (int y = 0; y < Ny; y++) {
for (int x = 0; x < Nx; x++) {
@ -776,12 +766,49 @@ void FlatSourceDomain::output_to_vtk() const
sample.x = ll.x + x_delta / 2.0 + x * x_delta;
Particle p;
p.r() = sample;
p.r_last() = sample;
p.E() = 1.0;
p.E_last() = 1.0;
p.u() = {1.0, 0.0, 0.0};
bool found = exhaustive_find_cell(p);
if (!found) {
voxel_indices[z * Ny * Nx + y * Nx + x] = -1;
voxel_positions[z * Ny * Nx + y * Nx + x] = sample;
weight_windows[z * Ny * Nx + y * Nx + x] = 0.0;
continue;
}
int i_cell = p.lowest_coord().cell;
int64_t source_region_idx =
source_region_offsets_[i_cell] + p.cell_instance();
voxel_indices[z * Ny * Nx + y * Nx + x] = source_region_idx;
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;
}
}
voxel_indices[z * Ny * Nx + y * Nx + x] = sr;
voxel_positions[z * Ny * Nx + y * Nx + x] = sample;
if (variance_reduction::weight_windows.size() == 1) {
WeightWindow ww =
variance_reduction::weight_windows[0]->get_weight_window(p);
float weight = ww.lower_weight;
weight_windows[z * Ny * Nx + y * Nx + x] = weight;
if (weight < min_weight)
min_weight = weight;
}
}
}
}
@ -798,10 +825,14 @@ void FlatSourceDomain::output_to_vtk() const
std::fprintf(plot, "BINARY\n");
std::fprintf(plot, "DATASET STRUCTURED_POINTS\n");
std::fprintf(plot, "DIMENSIONS %d %d %d\n", Nx, Ny, Nz);
std::fprintf(plot, "ORIGIN 0 0 0\n");
std::fprintf(plot, "ORIGIN %lf %lf %lf\n", ll.x, ll.y, ll.z);
std::fprintf(plot, "SPACING %lf %lf %lf\n", x_delta, y_delta, z_delta);
std::fprintf(plot, "POINT_DATA %d\n", Nx * Ny * Nz);
int64_t num_neg = 0;
int64_t num_samples = 0;
float min_flux = 0.0;
float max_flux = -1.0e20;
// Plot multigroup flux data
for (int g = 0; g < negroups_; g++) {
std::fprintf(plot, "SCALARS flux_group_%d float\n", g);
@ -809,12 +840,36 @@ void FlatSourceDomain::output_to_vtk() const
for (int i = 0; i < Nx * Ny * Nz; i++) {
int64_t fsr = voxel_indices[i];
int64_t source_element = fsr * negroups_ + g;
float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
float flux = 0;
if (fsr >= 0) {
flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
if (flux < 0.0)
flux = FlatSourceDomain::evaluate_flux_at_point(
voxel_positions[i], fsr, g);
}
if (flux < 0.0) {
num_neg++;
if (flux < min_flux) {
min_flux = flux;
}
}
if (flux > max_flux)
max_flux = flux;
num_samples++;
flux = convert_to_big_endian<float>(flux);
std::fwrite(&flux, sizeof(float), 1, plot);
}
}
// Slightly negative fluxes can be normal when sampling corners of linear
// source regions. However, very common and high magnitude negative fluxes
// may indicate numerical instability.
if (num_neg > 0) {
warning(fmt::format("{} plot samples ({:.4f}%) contained negative fluxes "
"(minumum found = {:.2e} maximum_found = {:.2e})",
num_neg, (100.0 * num_neg) / num_samples, min_flux, max_flux));
}
// Plot FSRs
std::fprintf(plot, "SCALARS FSRs float\n");
std::fprintf(plot, "LOOKUP_TABLE default\n");
@ -828,7 +883,9 @@ void FlatSourceDomain::output_to_vtk() const
std::fprintf(plot, "SCALARS Materials int\n");
std::fprintf(plot, "LOOKUP_TABLE default\n");
for (int fsr : voxel_indices) {
int mat = source_regions_.material(fsr);
int mat = -1;
if (fsr >= 0)
mat = source_regions_.material(fsr);
mat = convert_to_big_endian<int>(mat);
std::fwrite(&mat, sizeof(int), 1, plot);
}
@ -839,15 +896,16 @@ void FlatSourceDomain::output_to_vtk() const
std::fprintf(plot, "LOOKUP_TABLE default\n");
for (int i = 0; i < Nx * Ny * Nz; i++) {
int64_t fsr = voxel_indices[i];
float total_fission = 0.0;
int mat = source_regions_.material(fsr);
if (mat != MATERIAL_VOID) {
for (int g = 0; g < negroups_; g++) {
int64_t source_element = fsr * negroups_ + g;
float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
double sigma_f = sigma_f_[mat * negroups_ + g];
total_fission += sigma_f * flux;
if (fsr >= 0) {
int mat = source_regions_.material(fsr);
if (mat != MATERIAL_VOID) {
for (int g = 0; g < negroups_; g++) {
int64_t source_element = fsr * negroups_ + g;
float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g);
double sigma_f = sigma_f_[mat * negroups_ + g];
total_fission += sigma_f * flux;
}
}
}
total_fission = convert_to_big_endian<float>(total_fission);
@ -859,14 +917,29 @@ void FlatSourceDomain::output_to_vtk() const
for (int i = 0; i < Nx * Ny * Nz; i++) {
int64_t fsr = voxel_indices[i];
float total_external = 0.0f;
for (int g = 0; g < negroups_; g++) {
total_external += source_regions_.external_source(fsr, g);
if (fsr >= 0) {
for (int g = 0; g < negroups_; g++) {
total_external += source_regions_.external_source(fsr, g);
}
}
total_external = convert_to_big_endian<float>(total_external);
std::fwrite(&total_external, sizeof(float), 1, plot);
}
}
// Plot weight window data
if (variance_reduction::weight_windows.size() == 1) {
std::fprintf(plot, "SCALARS weight_window_lower float\n");
std::fprintf(plot, "LOOKUP_TABLE default\n");
for (int i = 0; i < Nx * Ny * Nz; i++) {
float weight = weight_windows[i];
if (weight == 0.0)
weight = min_weight;
weight = convert_to_big_endian<float>(weight);
std::fwrite(&weight, sizeof(float), 1, plot);
}
}
std::fclose(plot);
}
}
@ -938,7 +1011,7 @@ void FlatSourceDomain::count_external_source_regions()
{
n_external_source_regions_ = 0;
#pragma omp parallel for reduction(+ : n_external_source_regions_)
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
if (source_regions_.external_source_present(sr)) {
n_external_source_regions_++;
}
@ -984,7 +1057,7 @@ void FlatSourceDomain::convert_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++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
@ -1056,7 +1129,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
// small in magnitude, but rather due to the source region being physically
// small in volume and thus having a noisy flux estimate.
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
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) {
@ -1064,13 +1137,16 @@ void FlatSourceDomain::set_adjoint_sources(const vector<double>& forward_flux)
} else {
source_regions_.external_source(sr, g) = 1.0 / flux;
}
if (flux > 0.0) {
source_regions_.external_source_present(sr) = 1;
}
}
}
// 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++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
int material = source_regions_.material(sr);
if (material == MATERIAL_VOID) {
continue;
@ -1103,12 +1179,252 @@ void FlatSourceDomain::transpose_scattering_matrix()
void FlatSourceDomain::serialize_final_fluxes(vector<double>& flux)
{
// Ensure array is correct size
flux.resize(n_source_regions_ * negroups_);
flux.resize(n_source_regions() * negroups_);
// Serialize the final fluxes for output
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
flux[se] = source_regions_.scalar_flux_final(se);
}
}
void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell,
int32_t mesh_idx, int target_material_id, const vector<int32_t>& instances,
bool is_target_void)
{
Cell& cell = *model::cells[i_cell];
if (cell.type_ != Fill::MATERIAL)
return;
for (int32_t j : instances) {
int cell_material_idx = cell.material(j);
int cell_material_id = (cell_material_idx == C_NONE)
? C_NONE
: model::materials[cell_material_idx]->id();
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:
fatal_error(fmt::format("Source region {} already has mesh idx {} "
"applied, but trying to apply mesh idx {}",
sr, source_regions_.mesh(sr), mesh_idx));
}
source_regions_.mesh(sr) = mesh_idx;
}
}
}
void FlatSourceDomain::apply_mesh_to_cell_and_children(int32_t i_cell,
int32_t mesh_idx, int32_t target_material_id, bool is_target_void)
{
Cell& cell = *model::cells[i_cell];
if (cell.type_ == Fill::MATERIAL) {
vector<int> instances(cell.n_instances_);
std::iota(instances.begin(), instances.end(), 0);
apply_mesh_to_cell_instances(
i_cell, mesh_idx, target_material_id, instances, is_target_void);
} else if (target_material_id == C_NONE && !is_target_void) {
for (int j = 0; j < cell.n_instances_; j++) {
std::unordered_map<int32_t, vector<int32_t>> cell_instance_list =
cell.get_contained_cells(j, nullptr);
for (const auto& pair : cell_instance_list) {
int32_t i_child_cell = pair.first;
apply_mesh_to_cell_instances(i_child_cell, mesh_idx, target_material_id,
pair.second, is_target_void);
}
}
}
}
void FlatSourceDomain::apply_meshes()
{
// Skip if there are no mappings between mesh IDs and domains
if (mesh_domain_map_.empty())
return;
// Loop over meshes
for (int mesh_idx = 0; mesh_idx < model::meshes.size(); mesh_idx++) {
Mesh* mesh = model::meshes[mesh_idx].get();
int mesh_id = mesh->id();
// Skip if mesh id is not present in the map
if (mesh_domain_map_.find(mesh_id) == mesh_domain_map_.end())
continue;
// Loop over domains associated with the mesh
for (auto& domain : mesh_domain_map_[mesh_id]) {
Source::DomainType domain_type = domain.first;
int domain_id = domain.second;
if (domain_type == Source::DomainType::MATERIAL) {
for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) {
if (domain_id == C_NONE) {
apply_mesh_to_cell_and_children(i_cell, mesh_idx, domain_id, true);
} else {
apply_mesh_to_cell_and_children(i_cell, mesh_idx, domain_id, false);
}
}
} else if (domain_type == Source::DomainType::CELL) {
int32_t i_cell = model::cell_map[domain_id];
apply_mesh_to_cell_and_children(i_cell, mesh_idx, C_NONE, false);
} else if (domain_type == Source::DomainType::UNIVERSE) {
int32_t i_universe = model::universe_map[domain_id];
Universe& universe = *model::universes[i_universe];
for (int32_t i_cell : universe.cells_) {
apply_mesh_to_cell_and_children(i_cell, mesh_idx, C_NONE, false);
}
}
}
}
}
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 {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
// the source region key is found, we translate the key into a specific 1D
// source region index and return a handle its position in the
// source_regions_ vector.
auto it = source_region_map_.find(sr_key);
if (it != source_region_map_.end()) {
int64_t sr = it->second;
return source_regions_.get_source_region_handle(sr);
}
// Case 2: Check if the source region key is present in the temporary (thread
// safe) map. This is a common occurrence in the first power iteration when
// the source region has already been visited already by some other ray. We
// begin by locking the temporary map before any operations are performed. The
// lock is not global over the full data structure -- it will be dependent on
// which key is used.
discovered_source_regions_.lock(sr_key);
// If the key is found in the temporary map, then we return a handle to the
// source region that is stored in the temporary map.
if (discovered_source_regions_.contains(sr_key)) {
SourceRegionHandle handle {discovered_source_regions_[sr_key]};
discovered_source_regions_.unlock(sr_key);
return handle;
}
// Case 3: The source region key is not present anywhere, but it is only due
// to floating point artifacts. These artifacts occur when the overlaid mesh
// overlaps with actual geometry surfaces. In these cases, roundoff error may
// result in the ray tracer detecting an additional (very short) segment
// though a mesh bin that is actually past the physical source region
// boundary. This is a result of the the multi-level ray tracing treatment in
// OpenMC, which depending on the number of universes in the hierarchy etc can
// result in the wrong surface being selected as the nearest. This can happen
// in a lattice when there are two directions that both are very close in
// distance, within the tolerance of FP_REL_PRECISION, and the are thus
// treated as being equivalent so alternative logic is used. However, when we
// go and ray trace on this with the mesh tracer we may go past the surface
// bounding the current source region.
//
// To filter out this case, before we create the new source region, we double
// check that the actual starting point of this segment (r) is still in the
// same geometry source region that we started in. If an artifact is detected,
// we discard the segment (and attenuation through it) as it is not really a
// valid source region and will have only an infinitessimally small cell
// combined with the mesh bin. Thankfully, this is a fairly rare condition,
// and only triggers for very short ray lengths. It can be fixed by decreasing
// the value of FP_REL_PRECISION in constants.h, but this may have unknown
// consequences for the general ray tracer, so for now we do the below sanity
// checks before generating phantom source regions. A significant extra cost
// is incurred in instantiating the GeometryState object and doing a cell
// lookup, but again, this is going to be an extremely rare thing to check
// after the first power iteration has completed.
// Sanity check on source region id
GeometryState gs;
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) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
return handle;
}
// Sanity check on mesh bin
int mesh_idx = base_source_regions_.mesh(sr);
if (mesh_idx == C_NONE) {
if (mesh_bin != 0) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
return handle;
}
} else {
Mesh* mesh = model::meshes[mesh_idx].get();
int bin_found = mesh->get_bin(r + TINY_BIT * u);
if (bin_found != mesh_bin) {
discovered_source_regions_.unlock(sr_key);
SourceRegionHandle handle;
handle.is_numerical_fp_artifact_ = true;
return handle;
}
}
// Case 4: The source region key is valid, but is not present anywhere. This
// 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);
SourceRegionHandle handle {*sr_ptr};
return handle;
}
void FlatSourceDomain::finalize_discovered_source_regions()
{
// Extract keys for entries with a valid volume.
vector<SourceRegionKey> keys;
for (const auto& pair : discovered_source_regions_) {
if (pair.second.volume_ > 0.0) {
keys.push_back(pair.first);
}
}
if (!keys.empty()) {
// Sort the keys, so as to ensure reproducible ordering given that source
// regions may have been added to discovered_source_regions_ in an arbitrary
// order due to shared memory threading.
std::sort(keys.begin(), keys.end());
// Append the source regions in the sorted key order.
for (const auto& key : keys) {
const SourceRegion& sr = discovered_source_regions_[key];
source_region_map_[key] = source_regions_.n_source_regions();
source_regions_.push_back(sr);
}
// If any new source regions were discovered, we need to update the
// tally mapping between source regions and tally bins.
mapped_all_tallies_ = false;
}
discovered_source_regions_.clear();
}
} // namespace openmc

35
src/random_ray/linear_source_domain.cpp
View file

@ -24,12 +24,12 @@ void LinearSourceDomain::batch_reset()
{
FlatSourceDomain::batch_reset();
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
source_regions_.centroid_iteration(sr) = {0.0, 0.0, 0.0};
source_regions_.mom_matrix(sr) = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
}
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.flux_moments_new(se) = {0.0, 0.0, 0.0};
}
}
@ -40,8 +40,14 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
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 sr = 0; sr < n_source_regions_; sr++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) = 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;
@ -98,7 +104,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff)
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++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.source(se) += source_regions_.external_source(se);
}
}
@ -115,7 +121,7 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
// Normalize flux to total distance travelled by all rays this iteration
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.scalar_flux_new(se) *= normalization_factor;
source_regions_.flux_moments_new(se) *= normalization_factor;
}
@ -123,7 +129,7 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes(
// Accumulate cell-wise ray length tallies collected this iteration, then
// update the simulation-averaged cell-wise volume estimates
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr);
source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr);
source_regions_.volume_t(sr) += source_regions_.volume(sr);
@ -154,13 +160,22 @@ void LinearSourceDomain::set_flux_to_flux_plus_source(
source_regions_.scalar_flux_new(sr, g) /= volume;
source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
}
source_regions_.flux_moments_new(sr, g) *= (1.0 / volume);
// If a source region is small, then the moments are likely noisy, so we zero
// them. This is reasonable, given that small regions can get by with a flat
// source approximation anyhow.
if (source_regions_.is_small(sr)) {
source_regions_.flux_moments_new(sr, g) = {0.0, 0.0, 0.0};
} else {
source_regions_.flux_moments_new(sr, g) *= (1.0 / volume);
}
}
void LinearSourceDomain::set_flux_to_old_flux(int64_t sr, int g)
{
source_regions_.scalar_flux_new(g) = source_regions_.scalar_flux_old(g);
source_regions_.flux_moments_new(g) = source_regions_.flux_moments_old(g);
source_regions_.scalar_flux_new(sr, g) =
source_regions_.scalar_flux_old(sr, g);
source_regions_.flux_moments_new(sr, g) =
source_regions_.flux_moments_old(sr, g);
}
void LinearSourceDomain::accumulate_iteration_flux()
@ -170,7 +185,7 @@ void LinearSourceDomain::accumulate_iteration_flux()
// Accumulate scalar flux moments
#pragma omp parallel for
for (int64_t se = 0; se < n_source_elements_; se++) {
for (int64_t se = 0; se < n_source_elements(); se++) {
source_regions_.flux_moments_t(se) += source_regions_.flux_moments_new(se);
}
}

315
src/random_ray/random_ray.cpp
View file

@ -237,6 +237,7 @@ 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()
@ -313,7 +314,7 @@ void RandomRay::event_advance_ray()
wgt() = 0.0;
}
attenuate_flux(distance_alive, true);
attenuate_flux(distance_alive, true, distance_dead);
distance_travelled_ = distance_alive;
} else {
distance_travelled_ += distance;
@ -327,23 +328,84 @@ void RandomRay::event_advance_ray()
}
}
void RandomRay::attenuate_flux(double distance, bool is_active)
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();
// 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);
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();
}
}
} else {
attenuate_flux_inner(distance, is_active, sr, C_NONE, r());
}
}
void RandomRay::attenuate_flux_inner(
double distance, bool is_active, int64_t sr, int mesh_bin, Position r)
{
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);
}
switch (source_shape_) {
case RandomRaySourceShape::FLAT:
if (this->material() == MATERIAL_VOID) {
attenuate_flux_flat_source_void(distance, is_active);
attenuate_flux_flat_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_flat_source(distance, is_active);
attenuate_flux_flat_source(srh, distance, is_active, r);
}
break;
case RandomRaySourceShape::LINEAR:
case RandomRaySourceShape::LINEAR_XY:
if (this->material() == MATERIAL_VOID) {
attenuate_flux_linear_source_void(distance, is_active);
attenuate_flux_linear_source_void(srh, distance, is_active, r);
} else {
attenuate_flux_linear_source(distance, is_active);
attenuate_flux_linear_source(srh, distance, is_active, r);
}
break;
default:
@ -364,18 +426,13 @@ void RandomRay::attenuate_flux(double distance, bool is_active)
// than use of many atomic operations corresponding to each energy group
// individually (at least on CPU). Several other bookkeeping tasks are also
// performed when inside the lock.
void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
void RandomRay::attenuate_flux_flat_source(
SourceRegionHandle& srh, double distance, bool is_active, Position r)
{
// The number of geometric intersections is counted for reporting purposes
n_event()++;
// Determine source region index etc.
int i_cell = lowest_coord().cell;
// The source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
// The source element is the energy-specific region index
// Get material
int material = this->material();
// MOC incoming flux attenuation + source contribution/attenuation equation
@ -383,115 +440,99 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active)
float sigma_t = domain_->sigma_t_[material * negroups_ + g];
float tau = sigma_t * distance;
float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau)
float new_delta_psi =
(angular_flux_[g] - domain_->source_regions_.source(sr, g)) * exponential;
float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential;
delta_psi_[g] = new_delta_psi;
angular_flux_[g] -= new_delta_psi;
}
// If ray is in the active phase (not in dead zone), make contributions to
// source region bookkeeping
// Aquire lock for source region
srh.lock();
if (is_active) {
// Aquire lock for source region
domain_->source_regions_.lock(sr).lock();
// Accumulate delta psi into new estimate of source region flux for
// this iteration
for (int g = 0; g < negroups_; g++) {
domain_->source_regions_.scalar_flux_new(sr, g) += delta_psi_[g];
srh.scalar_flux_new(g) += delta_psi_[g];
}
// Accomulate volume (ray distance) into this iteration's estimate
// of the source region's volume
domain_->source_regions_.volume(sr) += distance;
srh.volume() += distance;
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!domain_->source_regions_.position_recorded(sr)) {
Position midpoint = r() + u() * (distance / 2.0);
domain_->source_regions_.position(sr) = midpoint;
domain_->source_regions_.position_recorded(sr) = 1;
}
// Release lock
domain_->source_regions_.lock(sr).unlock();
srh.n_hits() += 1;
}
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!srh.position_recorded()) {
Position midpoint = r + u() * (distance / 2.0);
srh.position() = midpoint;
srh.position_recorded() = 1;
}
// Release lock
srh.unlock();
}
// Alternative flux attenuation function for true void regions.
void RandomRay::attenuate_flux_flat_source_void(double distance, bool is_active)
void RandomRay::attenuate_flux_flat_source_void(
SourceRegionHandle& srh, double distance, bool is_active, Position r)
{
// The number of geometric intersections is counted for reporting purposes
n_event()++;
// Determine source region index etc.
int i_cell = lowest_coord().cell;
// The source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
int material = this->material();
// If ray is in the active phase (not in dead zone), make contributions to
// source region bookkeeping
if (is_active) {
// Aquire lock for source region
domain_->source_regions_.lock(sr).lock();
srh.lock();
// Accumulate delta psi into new estimate of source region flux for
// this iteration
for (int g = 0; g < negroups_; g++) {
domain_->source_regions_.scalar_flux_new(sr, g) +=
angular_flux_[g] * distance;
srh.scalar_flux_new(g) += angular_flux_[g] * distance;
}
// Accomulate volume (ray distance) into this iteration's estimate
// of the source region's volume
domain_->source_regions_.volume(sr) += distance;
domain_->source_regions_.volume_sq(sr) += distance * distance;
srh.volume() += distance;
srh.volume_sq() += distance * distance;
srh.n_hits() += 1;
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!domain_->source_regions_.position_recorded(sr)) {
Position midpoint = r() + u() * (distance / 2.0);
domain_->source_regions_.position(sr) = midpoint;
domain_->source_regions_.position_recorded(sr) = 1;
if (!srh.position_recorded()) {
Position midpoint = r + u() * (distance / 2.0);
srh.position() = midpoint;
srh.position_recorded() = 1;
}
// Release lock
domain_->source_regions_.lock(sr).unlock();
srh.unlock();
}
// Add source to incoming angular flux, assuming void region
for (int g = 0; g < negroups_; g++) {
angular_flux_[g] +=
domain_->source_regions_.external_source(sr, g) * distance;
angular_flux_[g] += srh.external_source(g) * distance;
}
}
void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
void RandomRay::attenuate_flux_linear_source(
SourceRegionHandle& srh, double distance, bool is_active, Position r)
{
// Cast domain to LinearSourceDomain
LinearSourceDomain* domain = dynamic_cast<LinearSourceDomain*>(domain_);
if (!domain) {
fatal_error("RandomRay::attenuate_flux_linear_source() called with "
"non-LinearSourceDomain domain.");
}
// The number of geometric intersections is counted for reporting purposes
n_event()++;
// Determine source region index etc.
int i_cell = lowest_coord().cell;
// The source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
// The source element is the energy-specific region index
int material = this->material();
Position& centroid = domain_->source_regions_.centroid(sr);
Position midpoint = r() + u() * (distance / 2.0);
Position& centroid = srh.centroid();
Position midpoint = r + u() * (distance / 2.0);
// Determine the local position of the midpoint and the ray origin
// relative to the source region's centroid
@ -503,9 +544,9 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
// be no estimate of its centroid. We detect this by checking if it has
// any accumulated volume. If its volume is zero, just use the midpoint
// of the ray as the region's centroid.
if (domain_->source_regions_.volume_t(sr)) {
if (srh.volume_t()) {
rm_local = midpoint - centroid;
r0_local = r() - centroid;
r0_local = r - centroid;
} else {
rm_local = {0.0, 0.0, 0.0};
r0_local = -u() * 0.5 * distance;
@ -530,10 +571,8 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
// calculated from the source gradients dot product with local centroid
// and direction, respectively.
float spatial_source =
domain_->source_regions_.source(sr, g) +
rm_local.dot(domain_->source_regions_.source_gradients(sr, g));
float dir_source =
u().dot(domain_->source_regions_.source_gradients(sr, g));
srh.source(g) + rm_local.dot(srh.source_gradients(g));
float dir_source = u().dot(srh.source_gradients(g));
float gn = exponentialG(tau);
float f1 = 1.0f - tau * gn;
@ -566,44 +605,47 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
}
}
// Compute an estimate of the spatial moments matrix for the source
// region based on parameters from this ray's crossing
MomentMatrix moment_matrix_estimate;
moment_matrix_estimate.compute_spatial_moments_matrix(
rm_local, u(), distance);
// Aquire lock for source region
srh.lock();
// If ray is in the active phase (not in dead zone), make contributions to
// source region bookkeeping
if (is_active) {
// Compute an estimate of the spatial moments matrix for the source
// region based on parameters from this ray's crossing
MomentMatrix moment_matrix_estimate;
moment_matrix_estimate.compute_spatial_moments_matrix(
rm_local, u(), distance);
// Aquire lock for source region
domain_->source_regions_.lock(sr).lock();
// Accumulate deltas into the new estimate of source region flux for this
// iteration
for (int g = 0; g < negroups_; g++) {
domain_->source_regions_.scalar_flux_new(sr, g) += delta_psi_[g];
domain_->source_regions_.flux_moments_new(sr, g) += delta_moments_[g];
srh.scalar_flux_new(g) += delta_psi_[g];
srh.flux_moments_new(g) += delta_moments_[g];
}
// Accumulate the volume (ray segment distance), centroid, and spatial
// momement estimates into the running totals for the iteration for this
// source region. The centroid and spatial momements estimates are scaled by
// the ray segment length as part of length averaging of the estimates.
domain_->source_regions_.volume(sr) += distance;
domain_->source_regions_.centroid_iteration(sr) += midpoint * distance;
// source region. The centroid and spatial momements estimates are scaled
// by the ray segment length as part of length averaging of the estimates.
srh.volume() += distance;
srh.centroid_iteration() += midpoint * distance;
moment_matrix_estimate *= distance;
domain_->source_regions_.mom_matrix(sr) += moment_matrix_estimate;
srh.mom_matrix() += moment_matrix_estimate;
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!domain_->source_regions_.position_recorded(sr)) {
domain_->source_regions_.position(sr) = midpoint;
domain_->source_regions_.position_recorded(sr) = 1;
}
// Release lock
domain_->source_regions_.lock(sr).unlock();
srh.n_hits() += 1;
}
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!srh.position_recorded()) {
srh.position() = midpoint;
srh.position_recorded() = 1;
}
// Release lock
srh.unlock();
}
// If traveling through a void region, the source term is either zero
@ -615,26 +657,13 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active)
// estimating the flux at specific pixel coordinates. Thus, plots will look
// nicer/more accurate if we record flux moments, so this function is useful.
void RandomRay::attenuate_flux_linear_source_void(
double distance, bool is_active)
SourceRegionHandle& srh, double distance, bool is_active, Position r)
{
// Cast domain to LinearSourceDomain
LinearSourceDomain* domain = dynamic_cast<LinearSourceDomain*>(domain_);
if (!domain) {
fatal_error("RandomRay::attenuate_flux_linear_source() called with "
"non-LinearSourceDomain domain.");
}
// The number of geometric intersections is counted for reporting purposes
n_event()++;
// Determine source region index etc.
int i_cell = lowest_coord().cell;
// The source region is the spatial region index
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
Position& centroid = domain_->source_regions_.centroid(sr);
Position midpoint = r() + u() * (distance / 2.0);
Position& centroid = srh.centroid();
Position midpoint = r + u() * (distance / 2.0);
// Determine the local position of the midpoint and the ray origin
// relative to the source region's centroid
@ -646,9 +675,9 @@ void RandomRay::attenuate_flux_linear_source_void(
// be no estimate of its centroid. We detect this by checking if it has
// any accumulated volume. If its volume is zero, just use the midpoint
// of the ray as the region's centroid.
if (domain_->source_regions_.volume_t(sr)) {
if (srh.volume_t()) {
rm_local = midpoint - centroid;
r0_local = r() - centroid;
r0_local = r - centroid;
} else {
rm_local = {0.0, 0.0, 0.0};
r0_local = -u() * 0.5 * distance;
@ -659,7 +688,7 @@ void RandomRay::attenuate_flux_linear_source_void(
// transport through a void region is greatly simplified. Here we
// compute the updated flux moments.
for (int g = 0; g < negroups_; g++) {
float spatial_source = domain_->source_regions_.source(sr, g);
float spatial_source = srh.external_source(g);
float new_delta_psi = (angular_flux_[g] - spatial_source) * distance;
float h1 = 0.5f;
h1 = h1 * angular_flux_[g];
@ -688,41 +717,41 @@ void RandomRay::attenuate_flux_linear_source_void(
rm_local, u(), distance);
// Aquire lock for source region
domain_->source_regions_.lock(sr).lock();
srh.lock();
// Accumulate delta psi into new estimate of source region flux for
// this iteration, and update flux momements
for (int g = 0; g < negroups_; g++) {
domain_->source_regions_.scalar_flux_new(sr, g) +=
angular_flux_[g] * distance;
domain_->source_regions_.flux_moments_new(sr, g) += delta_moments_[g];
srh.scalar_flux_new(g) += angular_flux_[g] * distance;
srh.flux_moments_new(g) += delta_moments_[g];
}
// Accumulate the volume (ray segment distance), centroid, and spatial
// momement estimates into the running totals for the iteration for this
// source region. The centroid and spatial momements estimates are scaled by
// the ray segment length as part of length averaging of the estimates.
domain_->source_regions_.volume(sr) += distance;
domain_->source_regions_.volume_sq(sr) += distance_2;
domain_->source_regions_.centroid_iteration(sr) += midpoint * distance;
srh.volume() += distance;
srh.volume_sq() += distance_2;
srh.centroid_iteration() += midpoint * distance;
moment_matrix_estimate *= distance;
domain_->source_regions_.mom_matrix(sr) += moment_matrix_estimate;
srh.mom_matrix() += moment_matrix_estimate;
// Tally valid position inside the source region (e.g., midpoint of
// the ray) if not done already
if (!domain_->source_regions_.position_recorded(sr)) {
domain_->source_regions_.position(sr) = midpoint;
domain_->source_regions_.position_recorded(sr) = 1;
if (!srh.position_recorded()) {
srh.position() = midpoint;
srh.position_recorded() = 1;
}
srh.n_hits() += 1;
// Release lock
domain_->source_regions_.lock(sr).unlock();
srh.unlock();
}
// Add source to incoming angular flux, assuming void region
for (int g = 0; g < negroups_; g++) {
angular_flux_[g] +=
domain_->source_regions_.external_source(sr, g) * distance;
angular_flux_[g] += srh.external_source(g) * distance;
}
}
@ -738,7 +767,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
wgt() = 1.0;
// set identifier for particle
id() = simulation::work_index[mpi::rank] + ray_id;
id() = ray_id;
// generate source site using sample method
SourceSite site;
@ -775,8 +804,26 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
int i_cell = lowest_coord().cell;
int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
for (int g = 0; g < negroups_; g++) {
angular_flux_[g] = domain_->source_regions_.source(sr, g);
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);
}
}
}

185
src/random_ray/random_ray_simulation.cpp
View file

@ -14,6 +14,7 @@
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/timer.h"
#include "openmc/weight_windows.h"
namespace openmc {
@ -48,13 +49,17 @@ void openmc_run_random_ray()
// 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 fixed sources, if present
sim.prepare_fixed_sources();
sim.apply_fixed_sources_and_mesh_domains();
// Begin main simulation timer
simulation::time_total.start();
@ -77,12 +82,13 @@ void openmc_run_random_ray()
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();
// Reduce variables across MPI ranks
sim.reduce_simulation_statistics();
// Output all simulation results
sim.output_simulation_results();
}
@ -107,7 +113,9 @@ void openmc_run_random_ray()
RandomRaySimulation adjoint_sim;
// Initialize adjoint fixed sources, if present
adjoint_sim.prepare_fixed_sources_adjoint(forward_flux);
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();
@ -127,9 +135,6 @@ void openmc_run_random_ray()
// Finalize OpenMC
openmc_simulation_finalize();
// Reduce variables across MPI ranks
adjoint_sim.reduce_simulation_statistics();
// Output all simulation results
adjoint_sim.output_simulation_results();
}
@ -321,12 +326,36 @@ void validate_random_ray_inputs()
#ifdef OPENMC_MPI
if (mpi::n_procs > 1) {
warning(
"Domain replication in random ray is supported, but suffers from poor "
"scaling of source all-reduce operations. Performance may severely "
"degrade beyond just a few MPI ranks. Domain decomposition may be "
"implemented in the future to provide efficient scaling.");
"MPI parallelism is not supported by the random ray solver. All work "
"will be performed by rank 0. Domain decomposition may be implemented in "
"the future to provide efficient MPI scaling.");
}
#endif
// Warn about instability resulting from linear sources in small regions
// 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 "
"generated by mesh subdivision. Negative sources may result in low "
"quality FW-CADIS weight windows. We recommend you use flat source mode "
"when generating weight windows with an overlaid mesh tally.");
}
}
void openmc_reset_random_ray()
{
FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
FlatSourceDomain::volume_normalized_flux_tallies_ = false;
FlatSourceDomain::adjoint_ = false;
FlatSourceDomain::mesh_domain_map_.clear();
RandomRay::ray_source_.reset();
RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
RandomRay::mesh_subdivision_enabled_ = false;
RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
}
//==============================================================================
@ -361,19 +390,29 @@ RandomRaySimulation::RandomRaySimulation()
domain_->flatten_xs();
}
void RandomRaySimulation::prepare_fixed_sources()
void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
{
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(
vector<double>& forward_flux)
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)
{
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);
}
}
@ -382,75 +421,93 @@ void RandomRaySimulation::simulate()
{
// Random ray power iteration loop
while (simulation::current_batch < settings::n_batches) {
// Initialize the current batch
initialize_batch();
initialize_generation();
// Reset total starting particle weight used for normalizing tallies
simulation::total_weight = 1.0;
// MPI not supported in random ray solver, so all work is done by rank 0
// TODO: Implement domain decomposition for MPI parallelism
if (mpi::master) {
// Update source term (scattering + fission)
domain_->update_neutron_source(k_eff_);
// Reset total starting particle weight used for normalizing tallies
simulation::total_weight = 1.0;
// Reset scalar fluxes, iteration volume tallies, and region hit flags to
// zero
domain_->batch_reset();
// Update source term (scattering + fission)
domain_->update_neutron_source(k_eff_);
// Start timer for transport
simulation::time_transport.start();
// 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
// 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();
// Transport sweep over all random rays for the iteration
#pragma omp parallel for schedule(dynamic) \
reduction(+ : total_geometric_intersections_)
for (int i = 0; i < simulation::work_per_rank; i++) {
RandomRay ray(i, domain_.get());
total_geometric_intersections_ +=
ray.transport_history_based_single_ray();
}
for (int i = 0; i < settings::n_particles; i++) {
RandomRay ray(i, domain_.get());
total_geometric_intersections_ +=
ray.transport_history_based_single_ray();
}
simulation::time_transport.stop();
simulation::time_transport.stop();
// If using multiple MPI ranks, perform all reduce on all transport results
domain_->all_reduce_replicated_source_regions();
// 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();
}
// Normalize scalar flux and update volumes
domain_->normalize_scalar_flux_and_volumes(
settings::n_particles * RandomRay::distance_active_);
// Normalize scalar flux and update volumes
domain_->normalize_scalar_flux_and_volumes(
settings::n_particles * RandomRay::distance_active_);
// Add source to scalar flux, compute number of FSR hits
int64_t n_hits = domain_->add_source_to_scalar_flux();
// Add source to scalar flux, compute number of FSR hits
int64_t n_hits = domain_->add_source_to_scalar_flux();
if (settings::run_mode == RunMode::EIGENVALUE) {
// Compute random ray k-eff
k_eff_ = domain_->compute_k_eff(k_eff_);
if (settings::run_mode == RunMode::EIGENVALUE) {
// Compute random ray k-eff
k_eff_ = domain_->compute_k_eff(k_eff_);
// Store random ray k-eff into OpenMC's native k-eff variable
global_tally_tracklength = k_eff_;
}
// Store random ray k-eff into OpenMC's native k-eff variable
global_tally_tracklength = k_eff_;
}
// Execute all tallying tasks, if this is an active batch
if (simulation::current_batch > settings::n_inactive) {
// Execute all tallying tasks, if this is an active batch
if (simulation::current_batch > settings::n_inactive) {
// Add this iteration's scalar flux estimate to final accumulated estimate
domain_->accumulate_iteration_flux();
// Add this iteration's scalar flux estimate to final accumulated
// estimate
domain_->accumulate_iteration_flux();
if (mpi::master) {
// Generate mapping between source regions and tallies
if (!domain_->mapped_all_tallies_) {
domain_->convert_source_regions_to_tallies();
}
// Use above mapping to contribute FSR flux data to appropriate tallies
// Use above mapping to contribute FSR flux data to appropriate
// tallies
domain_->random_ray_tally();
}
}
// Set phi_old = phi_new
domain_->flux_swap();
// Set phi_old = phi_new
domain_->flux_swap();
// Check for any obvious insabilities/nans/infs
instability_check(n_hits, k_eff_, avg_miss_rate_);
// Check for any obvious insabilities/nans/infs
instability_check(n_hits, k_eff_, avg_miss_rate_);
} // End MPI master work
// Finalize the current batch
finalize_generation();
@ -458,27 +515,13 @@ void RandomRaySimulation::simulate()
} // End random ray power iteration loop
}
void RandomRaySimulation::reduce_simulation_statistics()
{
// Reduce number of intersections
#ifdef OPENMC_MPI
if (mpi::n_procs > 1) {
uint64_t total_geometric_intersections_reduced = 0;
MPI_Reduce(&total_geometric_intersections_,
&total_geometric_intersections_reduced, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0,
mpi::intracomm);
total_geometric_intersections_ = total_geometric_intersections_reduced;
}
#endif
}
void RandomRaySimulation::output_simulation_results() const
{
// Print random ray results
if (mpi::master) {
print_results_random_ray(total_geometric_intersections_,
avg_miss_rate_ / settings::n_batches, negroups_,
domain_->n_source_regions_, domain_->n_external_source_regions_);
domain_->n_source_regions(), domain_->n_external_source_regions_);
if (model::plots.size() > 0) {
domain_->output_to_vtk();
}
@ -490,8 +533,8 @@ void RandomRaySimulation::output_simulation_results() const
void RandomRaySimulation::instability_check(
int64_t n_hits, double k_eff, double& avg_miss_rate) const
{
double percent_missed = ((domain_->n_source_regions_ - n_hits) /
static_cast<double>(domain_->n_source_regions_)) *
double percent_missed = ((domain_->n_source_regions() - n_hits) /
static_cast<double>(domain_->n_source_regions())) *
100.0;
avg_miss_rate += percent_missed;

229
src/random_ray/source_region.cpp
View file

@ -2,9 +2,36 @@
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/simulation.h"
namespace openmc {
//==============================================================================
// SourceRegionHandle implementation
//==============================================================================
SourceRegionHandle::SourceRegionHandle(SourceRegion& sr)
: negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_),
is_small_(&sr.is_small_), n_hits_(&sr.n_hits_),
is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_),
volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_),
volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_),
position_recorded_(&sr.position_recorded_),
external_source_present_(&sr.external_source_present_),
position_(&sr.position_), centroid_(&sr.centroid_),
centroid_iteration_(&sr.centroid_iteration_), centroid_t_(&sr.centroid_t_),
mom_matrix_(&sr.mom_matrix_), mom_matrix_t_(&sr.mom_matrix_t_),
volume_task_(&sr.volume_task_), mesh_(&sr.mesh_),
parent_sr_(&sr.parent_sr_), scalar_flux_old_(sr.scalar_flux_old_.data()),
scalar_flux_new_(sr.scalar_flux_new_.data()), source_(sr.source_.data()),
external_source_(sr.external_source_.data()),
scalar_flux_final_(sr.scalar_flux_final_.data()),
source_gradients_(sr.source_gradients_.data()),
flux_moments_old_(sr.flux_moments_old_.data()),
flux_moments_new_(sr.flux_moments_new_.data()),
flux_moments_t_(sr.flux_moments_t_.data()),
tally_task_(sr.tally_task_.data())
{}
//==============================================================================
// SourceRegion implementation
//==============================================================================
@ -25,7 +52,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear)
scalar_flux_final_.assign(negroups, 0.0);
tally_task_.resize(negroups);
if (is_linear) {
source_gradients_.resize(negroups);
flux_moments_old_.resize(negroups);
@ -34,15 +60,37 @@ 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
//==============================================================================
void SourceRegionContainer::push_back(const SourceRegion& sr)
{
n_source_regions_++;
// Scalar fields
material_.push_back(sr.material_);
is_small_.push_back(sr.is_small_);
n_hits_.push_back(sr.n_hits_);
lock_.push_back(sr.lock_);
volume_.push_back(sr.volume_);
volume_t_.push_back(sr.volume_t_);
@ -53,6 +101,8 @@ void SourceRegionContainer::push_back(const SourceRegion& sr)
external_source_present_.push_back(sr.external_source_present_);
position_.push_back(sr.position_);
volume_task_.push_back(sr.volume_task_);
mesh_.push_back(sr.mesh_);
parent_sr_.push_back(sr.parent_sr_);
// Only store these fields if is_linear_ is true
if (is_linear_) {
@ -92,6 +142,8 @@ void SourceRegionContainer::assign(
// Clear existing data
n_source_regions_ = 0;
material_.clear();
is_small_.clear();
n_hits_.clear();
lock_.clear();
volume_.clear();
volume_t_.clear();
@ -101,6 +153,8 @@ void SourceRegionContainer::assign(
position_recorded_.clear();
external_source_present_.clear();
position_.clear();
mesh_.clear();
parent_sr_.clear();
if (is_linear_) {
centroid_.clear();
@ -140,103 +194,88 @@ void SourceRegionContainer::flux_swap()
}
}
void SourceRegionContainer::mpi_sync_ranks(bool reduce_position)
SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr)
{
#ifdef OPENMC_MPI
SourceRegionHandle handle;
handle.negroups_ = negroups();
handle.material_ = &material(sr);
handle.is_small_ = &is_small(sr);
handle.n_hits_ = &n_hits(sr);
handle.is_linear_ = is_linear();
handle.lock_ = &lock(sr);
handle.volume_ = &volume(sr);
handle.volume_t_ = &volume_t(sr);
handle.volume_sq_ = &volume_sq(sr);
handle.volume_sq_t_ = &volume_sq_t(sr);
handle.volume_naive_ = &volume_naive(sr);
handle.position_recorded_ = &position_recorded(sr);
handle.external_source_present_ = &external_source_present(sr);
handle.position_ = &position(sr);
handle.volume_task_ = &volume_task(sr);
handle.mesh_ = &mesh(sr);
handle.parent_sr_ = &parent_sr(sr);
handle.scalar_flux_old_ = &scalar_flux_old(sr, 0);
handle.scalar_flux_new_ = &scalar_flux_new(sr, 0);
handle.source_ = &source(sr, 0);
handle.external_source_ = &external_source(sr, 0);
handle.scalar_flux_final_ = &scalar_flux_final(sr, 0);
handle.tally_task_ = &tally_task(sr, 0);
// The "position_recorded" variable needs to be allreduced (and maxed),
// as whether or not a cell was hit will affect some decisions in how the
// source is calculated in the next iteration so as to avoid dividing
// by zero. We take the max rather than the sum as the hit values are
// expected to be zero or 1.
MPI_Allreduce(MPI_IN_PLACE, position_recorded_.data(), 1, MPI_INT, MPI_MAX,
mpi::intracomm);
// The position variable is more complicated to reduce than the others,
// as we do not want the sum of all positions in each cell, rather, we
// want to just pick any single valid position. Thus, we perform a gather
// and then pick the first valid position we find for all source regions
// that have had a position recorded. This operation does not need to
// be broadcast back to other ranks, as this value is only used for the
// tally conversion operation, which is only performed on the master rank.
// While this is expensive, it only needs to be done for active batches,
// and only if we have not mapped all the tallies yet. Once tallies are
// fully mapped, then the position vector is fully populated, so this
// operation can be skipped.
// Then, we perform the gather of position data, if needed
if (reduce_position) {
// Master rank will gather results and pick valid positions
if (mpi::master) {
// Initialize temporary vector for receiving positions
vector<vector<Position>> all_position;
all_position.resize(mpi::n_procs);
for (int i = 0; i < mpi::n_procs; i++) {
all_position[i].resize(n_source_regions_);
}
// Copy master rank data into gathered vector for convenience
all_position[0] = position_;
// Receive all data into gather vector
for (int i = 1; i < mpi::n_procs; i++) {
MPI_Recv(all_position[i].data(), n_source_regions_ * 3, MPI_DOUBLE, i,
0, mpi::intracomm, MPI_STATUS_IGNORE);
}
// Scan through gathered data and pick first valid cell posiiton
for (int64_t sr = 0; sr < n_source_regions_; sr++) {
if (position_recorded_[sr] == 1) {
for (int i = 0; i < mpi::n_procs; i++) {
if (all_position[i][sr].x != 0.0 || all_position[i][sr].y != 0.0 ||
all_position[i][sr].z != 0.0) {
position_[sr] = all_position[i][sr];
break;
}
}
}
}
} else {
// Other ranks just send in their data
MPI_Send(position_.data(), n_source_regions_ * 3, MPI_DOUBLE, 0, 0,
mpi::intracomm);
}
if (handle.is_linear_) {
handle.centroid_ = &centroid(sr);
handle.centroid_iteration_ = &centroid_iteration(sr);
handle.centroid_t_ = &centroid_t(sr);
handle.mom_matrix_ = &mom_matrix(sr);
handle.mom_matrix_t_ = &mom_matrix_t(sr);
handle.source_gradients_ = &source_gradients(sr, 0);
handle.flux_moments_old_ = &flux_moments_old(sr, 0);
handle.flux_moments_new_ = &flux_moments_new(sr, 0);
handle.flux_moments_t_ = &flux_moments_t(sr, 0);
}
// For the rest of the source region data, we simply perform an all reduce,
// as these values will be needed on all ranks for transport during the
// next iteration.
MPI_Allreduce(MPI_IN_PLACE, volume_.data(), n_source_regions_, MPI_DOUBLE,
MPI_SUM, mpi::intracomm);
MPI_Allreduce(MPI_IN_PLACE, volume_sq_.data(), n_source_regions_, MPI_DOUBLE,
MPI_SUM, mpi::intracomm);
MPI_Allreduce(MPI_IN_PLACE, scalar_flux_new_.data(),
n_source_regions_ * negroups_, MPI_DOUBLE, MPI_SUM, mpi::intracomm);
if (is_linear_) {
// We are going to assume we can safely cast Position, MomentArray,
// and MomentMatrix to contiguous arrays of doubles for the MPI
// allreduce operation. This is a safe assumption as typically
// compilers will at most pad to 8 byte boundaries. If a new FP32
// MomentArray type is introduced, then there will likely be padding, in
// which case this function will need to become more complex.
if (sizeof(MomentArray) != 3 * sizeof(double) ||
sizeof(MomentMatrix) != 6 * sizeof(double)) {
fatal_error(
"Unexpected buffer padding in linear source domain reduction.");
}
MPI_Allreduce(MPI_IN_PLACE, static_cast<void*>(flux_moments_new_.data()),
n_source_regions_ * negroups_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm);
MPI_Allreduce(MPI_IN_PLACE, static_cast<void*>(mom_matrix_.data()),
n_source_regions_ * 6, MPI_DOUBLE, MPI_SUM, mpi::intracomm);
MPI_Allreduce(MPI_IN_PLACE, static_cast<void*>(centroid_iteration_.data()),
n_source_regions_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm);
}
#endif
return handle;
}
} // namespace openmc
void SourceRegionContainer::adjoint_reset()
{
std::fill(n_hits_.begin(), n_hits_.end(), 0);
std::fill(volume_.begin(), volume_.end(), 0.0);
std::fill(volume_t_.begin(), volume_t_.end(), 0.0);
std::fill(volume_sq_.begin(), volume_sq_.end(), 0.0);
std::fill(volume_sq_t_.begin(), volume_sq_t_.end(), 0.0);
std::fill(volume_naive_.begin(), volume_naive_.end(), 0.0);
std::fill(
external_source_present_.begin(), external_source_present_.end(), 0);
std::fill(external_source_.begin(), external_source_.end(), 0.0);
std::fill(centroid_.begin(), centroid_.end(), Position {0.0, 0.0, 0.0});
std::fill(centroid_iteration_.begin(), centroid_iteration_.end(),
Position {0.0, 0.0, 0.0});
std::fill(centroid_t_.begin(), centroid_t_.end(), Position {0.0, 0.0, 0.0});
std::fill(mom_matrix_.begin(), mom_matrix_.end(),
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});
for (auto& task_set : volume_task_) {
task_set.clear();
}
std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.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(),
MomentArray {0.0, 0.0, 0.0});
std::fill(flux_moments_old_.begin(), flux_moments_old_.end(),
MomentArray {0.0, 0.0, 0.0});
std::fill(flux_moments_new_.begin(), flux_moments_new_.end(),
MomentArray {0.0, 0.0, 0.0});
std::fill(flux_moments_t_.begin(), flux_moments_t_.end(),
MomentArray {0.0, 0.0, 0.0});
for (auto& task_set : tally_task_) {
task_set.clear();
}
}
} // namespace openmc

26
src/settings.cpp
View file

@ -1,4 +1,5 @@
#include "openmc/settings.h"
#include "openmc/random_ray/flat_source_domain.h"
#include <cmath> // for ceil, pow
#include <limits> // for numeric_limits
@ -319,6 +320,31 @@ void get_run_parameters(pugi::xml_node node_base)
fatal_error("Unrecognized sample method: " + temp_str);
}
}
if (check_for_node(random_ray_node, "source_region_meshes")) {
pugi::xml_node node_source_region_meshes =
random_ray_node.child("source_region_meshes");
for (pugi::xml_node node_mesh :
node_source_region_meshes.children("mesh")) {
int mesh_id = std::stoi(node_mesh.attribute("id").value());
for (pugi::xml_node node_domain : node_mesh.children("domain")) {
int domain_id = std::stoi(node_domain.attribute("id").value());
std::string domain_type = node_domain.attribute("type").value();
Source::DomainType type;
if (domain_type == "material") {
type = Source::DomainType::MATERIAL;
} else if (domain_type == "cell") {
type = Source::DomainType::CELL;
} else if (domain_type == "universe") {
type = Source::DomainType::UNIVERSE;
} else {
throw std::runtime_error("Unknown domain type: " + domain_type);
}
FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back(
type, domain_id);
RandomRay::mesh_subdivision_enabled_ = true;
}
}
}
}
}

13
src/weight_windows.cpp
View file

@ -657,11 +657,18 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
}
}
xt::noalias(new_bounds) = 1.0 / new_bounds;
// We take the inverse, but are careful not to divide by zero e.g. if some
// mesh bins are not reachable in the physical geometry.
xt::noalias(new_bounds) =
xt::where(xt::not_equal(new_bounds, 0.0), 1.0 / new_bounds, 0.0);
auto max_val = xt::amax(new_bounds)();
xt::noalias(new_bounds) = new_bounds / (2.0 * max_val);
// For bins that were missed, we use the minimum weight window value. This
// shouldn't matter except for plotting.
auto min_val = xt::amin(new_bounds)();
xt::noalias(new_bounds) =
xt::where(xt::not_equal(new_bounds, 0.0), new_bounds, min_val);
}
// make sure that values where the mean is zero are set s.t. the weight window

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

271
tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat Normal file
View file

@ -0,0 +1,271 @@
<?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">
<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 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="1" type="universe"/>
</mesh>
<mesh id="2">
<domain id="2" type="cell"/>
</mesh>
<mesh id="3">
<domain id="3" type="material"/>
</mesh>
</source_region_meshes>
<source_shape>flat</source_shape>
</random_ray>
<mesh id="1">
<dimension>24 24 24</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<mesh id="2">
<dimension>36 36 36</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<mesh id="3">
<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>

9
tests/regression_tests/random_ray_fixed_source_mesh/flat/results_true.dat Normal file
View file

@ -0,0 +1,9 @@
tally 1:
1.750296E+00
2.051829E-01
tally 2:
8.045199E-02
4.384408E-04
tally 3:
5.216828E-03
1.840861E-06

271
tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat Normal file
View file

@ -0,0 +1,271 @@
<?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">
<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 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="1" type="universe"/>
</mesh>
<mesh id="2">
<domain id="2" type="cell"/>
</mesh>
<mesh id="3">
<domain id="3" type="material"/>
</mesh>
</source_region_meshes>
<source_shape>linear</source_shape>
</random_ray>
<mesh id="1">
<dimension>24 24 24</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<mesh id="2">
<dimension>36 36 36</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<mesh id="3">
<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>

9
tests/regression_tests/random_ray_fixed_source_mesh/linear/results_true.dat Normal file
View file

@ -0,0 +1,9 @@
tally 1:
1.780361E+00
2.137912E-01
tally 2:
8.230400E-02
4.596391E-04
tally 3:
5.207797E-03
1.834531E-06

53
tests/regression_tests/random_ray_fixed_source_mesh/test.py Normal file
View file

@ -0,0 +1,53 @@
import os
import openmc
from openmc.examples import random_ray_three_region_cube
from openmc.utility_funcs import change_directory
import pytest
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 make_mesh(dim):
width = 30.0
mesh = openmc.RegularMesh()
mesh.dimension = (dim, dim, dim)
mesh.lower_left = (0.0, 0.0, 0.0)
mesh.upper_right = (width, width, width)
return mesh
@pytest.mark.parametrize("shape", ["flat", "linear"])
def test_random_ray_fixed_source_mesh(shape):
with change_directory(shape):
openmc.reset_auto_ids()
model = random_ray_three_region_cube()
# We will apply three different mesh resolutions to three different domain types
source_universe = model.geometry.get_universes_by_name('source universe')[0]
void_cell = model.geometry.get_cells_by_name('infinite void region')[0]
absorber_mat = model.geometry.get_materials_by_name('absorber')[0]
model.settings.random_ray['source_region_meshes'] = [
(make_mesh(24), [source_universe]),
(make_mesh(36), [void_cell]),
(make_mesh(30), [absorber_mat])
]
# We also test flat/linear source shapes to ensure they are both
# working correctly with the mesh overlay logic
model.settings.random_ray['source_shape'] = shape
model.settings.inactive = 15
model.settings.batches = 30
harness = MGXSTestHarness('statepoint.30.h5', model)
harness.main()

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

119
tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat Normal file
View file

@ -0,0 +1,119 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0"/>
<macroscopic name="UO2"/>
</material>
<material id="2" name="Water">
<density units="macro" value="1.0"/>
<macroscopic name="LWTR"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
<cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
<cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
<cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
<cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
<cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
<cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
<cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
<cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
<cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
<cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
<cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
<cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
<cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
<cell id="15" material="2" name="moderator infinite" universe="3"/>
<cell fill="4" id="16" universe="5"/>
<cell fill="6" id="17" name="assembly" region="14 -15 16 -17" universe="7"/>
<lattice id="4">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</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 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 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>
<lattice id="6">
<pitch>1.26 1.26</pitch>
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<universes>
2 2
2 5 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
<surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
<surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
<surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
<surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
<surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
<surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
<surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
<surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
<surface boundary="reflective" coeffs="-1.26" id="14" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="1.26" id="15" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.26" id="16" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="1.26" id="17" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>100.0</distance_active>
<distance_inactive>20.0</distance_inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>True</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="2">
<domain id="7" type="universe"/>
</mesh>
</source_region_meshes>
</random_ray>
<mesh id="2">
<dimension>40 40</dimension>
<lower_left>-1.26 -1.26</lower_left>
<upper_right>1.26 1.26</upper_right>
</mesh>
</settings>
<tallies>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<upper_right>1.26 1.26</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<tally id="1" name="Mesh tally">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

171
tests/regression_tests/random_ray_k_eff_mesh/results_true.dat Normal file
View file

@ -0,0 +1,171 @@
k-combined:
8.379203E-01 8.057199E-03
tally 1:
5.080171E+00
5.167984E+00
1.880341E+00
7.079266E-01
4.576373E+00
4.193319E+00
2.859914E+00
1.638769E+00
4.243332E-01
3.607732E-02
1.032742E+00
2.136998E-01
1.692643E+00
5.794069E-01
5.445214E-02
5.995212E-04
1.325256E-01
3.551193E-03
2.372336E+00
1.147031E+00
7.807378E-02
1.242019E-03
1.900160E-01
7.356955E-03
7.135636E+00
1.035026E+01
8.273225E-02
1.392069E-03
2.013562E-01
8.245961E-03
2.044034E+01
8.394042E+01
3.100485E-02
1.932097E-04
7.671933E-02
1.182985E-03
1.313652E+01
3.451847E+01
1.764978E-01
6.230420E-03
4.909196E-01
4.820140E-02
7.585874E+00
1.150936E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.386790E+00
2.295327E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.820058E+00
6.729238E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.694013E+00
1.481363E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.453818E+00
1.128202E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.823642E+01
6.682239E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.137903E+01
2.590265E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.589529E+00
4.219985E+00
1.712446E+00
5.873675E-01
4.167750E+00
3.479202E+00
2.728285E+00
1.492132E+00
4.104063E-01
3.377132E-02
9.988468E-01
2.000404E-01
1.660506E+00
5.567967E-01
5.427391E-02
5.944708E-04
1.320918E-01
3.521277E-03
2.305889E+00
1.082691E+00
7.695793E-02
1.205644E-03
1.873002E-01
7.141490E-03
7.076637E+00
1.017946E+01
8.323139E-02
1.408687E-03
2.025710E-01
8.344398E-03
2.095897E+01
8.825741E+01
3.236513E-02
2.104220E-04
8.008525E-02
1.288373E-03
1.358006E+01
3.689205E+01
1.862562E-01
6.941428E-03
5.180621E-01
5.370209E-02
5.067386E+00
5.141748E+00
1.912704E+00
7.328484E-01
4.655140E+00
4.340941E+00
2.858992E+00
1.637705E+00
4.331050E-01
3.759875E-02
1.054091E+00
2.227118E-01
1.692974E+00
5.796396E-01
5.560487E-02
6.246120E-04
1.353311E-01
3.699815E-03
2.368737E+00
1.143184E+00
7.950085E-02
1.286135E-03
1.934892E-01
7.618268E-03
7.119767E+00
1.030095E+01
8.427627E-02
1.442764E-03
2.051141E-01
8.546249E-03
2.047651E+01
8.426275E+01
3.183786E-02
2.037156E-04
7.878057E-02
1.247311E-03
1.326415E+01
3.519010E+01
1.833688E-01
6.726832E-03
5.100312E-01
5.204188E-02

36
tests/regression_tests/random_ray_k_eff_mesh/test.py Normal file
View file

@ -0,0 +1,36 @@
import os
import openmc
from openmc.examples import random_ray_lattice
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_k_eff_mesh():
model = random_ray_lattice()
# The model already has some geometrical subdivisions
# up to a 10x10 grid in the moderator region. So, we
# increase the resolution 40x40 applied over the full
# 2x2 lattice.
pitch = 1.26
dim = 40
mesh = openmc.RegularMesh()
mesh.dimension = (dim, dim)
mesh.lower_left = (-pitch, -pitch)
mesh.upper_right = (pitch, pitch)
root = model.geometry.root_universe
model.settings.random_ray['source_region_meshes'] = [(mesh, [root])]
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()

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

265
tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat Normal file
View file

@ -0,0 +1,265 @@
<?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>750</particles>
<batches>30</batches>
<inactive>20</inactive>
<source particle="neutron" strength="3.14" type="independent">
<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>
<weight_window_generators>
<weight_windows_generator>
<mesh>1</mesh>
<particle_type>neutron</particle_type>
<max_realizations>10</max_realizations>
<update_interval>1</update_interval>
<on_the_fly>true</on_the_fly>
<method>fw_cadis</method>
</weight_windows_generator>
</weight_window_generators>
<mesh id="1">
<dimension>15 15 15</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<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>
<source_shape>flat</source_shape>
</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>

6760
tests/regression_tests/weightwindows_fw_cadis_mesh/flat/results_true.dat Normal file
View file

File diff suppressed because it is too large Load diff

265
tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat Normal file
View file

@ -0,0 +1,265 @@
<?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>750</particles>
<batches>30</batches>
<inactive>20</inactive>
<source particle="neutron" strength="3.14" type="independent">
<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>
<weight_window_generators>
<weight_windows_generator>
<mesh>1</mesh>
<particle_type>neutron</particle_type>
<max_realizations>10</max_realizations>
<update_interval>1</update_interval>
<on_the_fly>true</on_the_fly>
<method>fw_cadis</method>
</weight_windows_generator>
</weight_window_generators>
<mesh id="1">
<dimension>15 15 15</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>30.0 30.0 30.0</upper_right>
</mesh>
<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>
<source_shape>linear</source_shape>
</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>

6760
tests/regression_tests/weightwindows_fw_cadis_mesh/linear/results_true.dat Normal file
View file

File diff suppressed because it is too large Load diff

49
tests/regression_tests/weightwindows_fw_cadis_mesh/test.py Normal file
View file

@ -0,0 +1,49 @@
import os
import openmc
from openmc.utility_funcs import change_directory
from openmc.examples import random_ray_three_region_cube
import pytest
from tests.testing_harness import WeightWindowPyAPITestHarness
class MGXSTestHarness(WeightWindowPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
@pytest.mark.parametrize("shape", ["flat", "linear"])
def test_weight_windows_fw_cadis_mesh(shape):
with change_directory(shape):
openmc.reset_auto_ids()
model = random_ray_three_region_cube()
# The base model has a resolution of 12, so we overlay
# something else for FW-CADIS
n = 15
width = 30.0
ww_mesh = openmc.RegularMesh()
ww_mesh.dimension = (n, n, n)
ww_mesh.lower_left = (0.0, 0.0, 0.0)
ww_mesh.upper_right = (width, width, width)
wwg = openmc.WeightWindowGenerator(
method="fw_cadis", mesh=ww_mesh, max_realizations=model.settings.batches)
model.settings.weight_window_generators = wwg
root = model.geometry.root_universe
model.settings.random_ray['source_region_meshes'] = [(ww_mesh, [root])]
model.settings.particles = 750
model.settings.batches = 30
model.settings.inactive = 20
model.settings.random_ray['source_shape'] = shape
harness = MGXSTestHarness('statepoint.30.h5', model)
harness.main()