diff --git a/CMakeLists.txt b/CMakeLists.txt index 3cbeaa0241b..0f4cc1b527b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -380,6 +380,8 @@ list(APPEND libopenmc_SOURCES src/random_ray/random_ray_simulation.cpp src/random_ray/random_ray.cpp src/random_ray/flat_source_domain.cpp + src/random_ray/linear_source_domain.cpp + src/random_ray/moment_matrix.cpp src/reaction.cpp src/reaction_product.cpp src/scattdata.cpp diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 45af2b0c894..9f8eb84d80e 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -218,9 +218,9 @@ Following the multigroup discretization, another assumption made is that a large and complex problem can be broken up into small constant cross section regions, and that these regions have group dependent, flat, isotropic sources (fission and scattering), :math:`Q_g`. Anisotropic as well as higher order sources are -also possible with MOC-based methods but are not used yet in OpenMC for -simplicity. With these key assumptions, the multigroup MOC form of the neutron -transport equation can be written as in Equation :eq:`moc_final`. +also possible with MOC-based methods. With these key assumptions, the multigroup +MOC form of the neutron transport equation can be written as in Equation +:eq:`moc_final`. .. math:: :label: moc_final @@ -287,7 +287,7 @@ final expression for the average angular flux for a ray crossing a region as: .. math:: :label: average_psi_final - \overline{\psi}_{r,i,g} = \frac{Q_{i,g}}{\Sigma_{t,i,g}} + \frac{\Delta \psi_{r,g}}{\ell_r \Sigma_{t,i,g}} + \overline{\psi}_{r,i,g} = \frac{Q_{i,g}}{\Sigma_{t,i,g}} + \frac{\Delta \psi_{r,g}}{\ell_r \Sigma_{t,i,g}}. ~~~~~~~~~~~ Random Rays @@ -771,6 +771,170 @@ By default, the unnormalized flux values (units of cm) will be reported. If the user wishes to received volume normalized flux tallies, then an option for this is available, as described in the :ref:`User Guide`. +-------------- +Linear Sources +-------------- + +Instead of making a flat source approximation, as in the previous section, a +Linear Source (LS) approximation can be used. Different LS approximations have +been developed; the OpenMC implementation follows the MOC LS scheme described by +`Ferrer `_. The LS source along a characteristic is given by: + +.. math:: + :label: linear_source + + Q_{i,g}(s) = \bar{Q}_{r,i,g} + \hat{Q}_{r,i,g}(s-\ell_{r}/2), + +where the source, :math:`Q_{i,g}(s)`, varies linearly along the track and +:math:`\bar{Q}_{r,i,g}` and :math:`\hat{Q}_{r,i,g}` are track specific source +terms to define shortly. Integrating the source, as done in Equation +:eq:`moc_final`, leads to + +.. math:: + :label: lsr_attenuation + + \psi^{out}_{r,g}=\psi^{in}_{r,g} + \left(\frac{\bar{Q}_{r, i, g}}{\Sigma_{\mathrm{t}, i, g}}-\psi^{in}_{r,g}\right) + F_{1}\left(\tau_{i,g}\right)+\frac{\hat{Q}_{r, i, g}^{g}}{2\left(\Sigma_{\mathrm{t}, i,g}\right)^{2}} F_{2}\left(\tau_{i,g}\right), + +where for simplicity the term :math:`\tau_{i,g}` and the expoentials :math:`F_1` +and :math:`F_2` are introduced, given by: + +.. math:: + :label: tau + + \tau_{i,g} = \Sigma_{\mathrm{t,i,g}} \ell_{r} + +.. math:: + :label: f1 + + F_1(\tau) = 1 - e^{-\tau}, + +and + +.. math:: + :label: f2 + + F_{2}\left(\tau\right) = 2\left[\tau-F_{1}\left(\tau\right)\right]-\tau F_{1}\left(\tau\right). + + +To solve for the track specific source terms in Equation :eq:`linear_source` we +first define a local reference frame. If we now refer to :math:`\mathbf{r}` as +the global coordinate and introduce the source region specific coordinate +:math:`\mathbf{u}` such that, + +.. math:: + :label: local_coord + + \mathbf{u}_{r} = \mathbf{r}-\mathbf{r}_{\mathrm{c}}, + +where :math:`\mathbf{r}_{\mathrm{c}}` is the centroid of the source region of +interest. In turn :math:`\mathbf{u}_{r,\mathrm{c}}` and :math:`\mathbf{u}_{r,0}` +are the local centroid and entry positions of a ray. The computation of the +local and global centroids are described further by `Gunow `_. + +Using the local position, the source in a source region is given by: + +.. math:: + :label: region_source + + \tilde{Q}(\boldsymbol{x}) ={Q}_{i,g}+ \boldsymbol{\vec{Q}}_{i,g} \cdot \mathbf{u}_{r}\;\mathrm{,} + +This definition allows us to solve for our characteric source terms resulting in: + +.. math:: + :label: source_term_1 + + \bar{Q}_{r, i, g} = Q_{i,g} + \left[\mathbf{u}_{r,\mathrm{c}} \cdot \boldsymbol{\vec{Q}}_{i,g}\right], + +.. math:: + :label: source_term_2 + + \hat{Q}_{r, i, g} = \left[\boldsymbol{\Omega} \cdot \boldsymbol{\vec{Q}}_{i,g}\right]\;\mathrm{,} + +:math:`\boldsymbol{\Omega}` being the direction vector of the ray. The next step +is to solve for the LS source vector :math:`\boldsymbol{\vec{Q}}_{i,g}`. A +relationship between the LS source vector and the source moments, +:math:`\boldsymbol{\vec{q}}_{i,g}` can be derived, as in `Ferrer +`_ and `Gunow `_: + +.. math:: + :label: m_equation + + \mathbf{M}_{i} \boldsymbol{\vec{Q}}_{i,g} = \boldsymbol{\vec{q}}_{i,g} \;\mathrm{.} + +The spatial moments matrix :math:`M_i` in region :math:`i` represents the +spatial distribution of the 3D object composing the `source region +`_. This matrix is independent of the material of the source +region, fluxes, and any transport effects -- it is a purely geometric quantity. +It is a symmetric :math:`3\times3` matrix. While :math:`M_i` is not known +apriori to the simulation, similar to the source region volume, it can be +computed "on-the-fly" as a byproduct of the random ray integration process. Each +time a ray randomly crosses the region within its active length, an estimate of +the spatial moments matrix can be computed by using the midpoint of the ray as +an estimate of the centroid, and the distance and direction of the ray can be +used to inform the other spatial moments within the matrix. As this information +is purely geometric, the stochastic estimate of the centroid and spatial moments +matrix can be accumulated and improved over the entire duration of the +simulation, converging towards their true quantities. + +With an estimate of the spatial moments matrix :math:`M_i` resulting from the +ray tracing process naturally, the LS source vector +:math:`\boldsymbol{\vec{Q}}_{i,g}` can be obtained via a linear solve of +:eq:`m_equation`, or by the direct inversion of :math:`M_i`. However, to +accomplish this, we must first know the source moments +:math:`\boldsymbol{\vec{q}}_{i,g}`. Fortunately, the source moments are also +defined by the definition of the source: + +.. math:: + :label: source_moments + + q_{v, i, g}= \frac{\chi_{i,g}}{k_{eff}} \sum_{g^{\prime}=1}^{G} \nu + \Sigma_{\mathrm{f},i, g^{\prime}} \hat{\phi}_{v, i, g^{\prime}} + \sum_{g^{\prime}=1}^{G} + \Sigma_{\mathrm{s}, i, g^{\prime}\rightarrow g} \hat{\phi}_{v, i, g^{\prime}}\quad \forall v \in(x, y, z)\;\mathrm{,} + +where :math:`v` indicates the direction vector component, and we have introduced +the scalar flux moments :math:`\hat{\phi}`. The scalar flux moments can be +solved for by taking the `integral definition `_ of a spatial +moment, allowing us to derive a "simulation averaged" estimator for the scalar +moment, as in Equation :eq:`phi_sim`, + +.. math:: + :label: scalar_moments_sim + + \hat{\phi}_{v,i,g}^{simulation} = \frac{\sum\limits_{r=1}^{N_i} + \ell_{r} \left[\Omega_{v} \hat{\psi}_{r,i,g} + u_{r,v,0} \bar{\psi}_{r,i,g}\right]} + {\Sigma_{t,i,g} \frac{\sum\limits^{B}_{b}\sum\limits^{N_i}_{r} \ell_{b,r} }{B}} + \quad \forall v \in(x, y, z)\;\mathrm{,} + + +where the average angular flux is given by Equation :eq:`average_psi_final`, and +the angular flux spatial moments :math:`\hat{\psi}_{r,i,g}` by: + +.. math:: + :label: angular_moments + + \hat{\psi}_{r, i, g} = \frac{\ell_{r}\psi^{in}_{r,g}}{2} + + \left(\frac{\bar{Q}_{r,i, g}}{\Sigma_{\mathrm{t}, i, g}}-\psi^{in}_{r,g}\right) + \frac{G_{1}\left(\tau_{i,g}\right)}{\Sigma_{\mathrm{t}, i, g}} + \frac{\ell_{r}\hat{Q}_{r,i,g}} + {2\left(\Sigma_{\mathrm{t}, i, g}\right)^{2}}G_{2}\left(\tau_{i,g}\right)\;\mathrm{.} + + +The new exponentials introduced, again for simplicity, are simply: + +.. math:: + :label: G1 + + G_{1}(\tau) = 1+\frac{\tau}{2}-\left(1+\frac{1}{\tau}\right) F_{1}(\tau), + +.. math:: + :label: G2 + + G_{2}(\tau) = \frac{2}{3} \tau-\left(1+\frac{2}{\tau}\right) G_{1}(\tau) + +The contents of this section, alongside the equations for the flat source and +scalar flux, Equations :eq:`source_update` and :eq:`phi_sim` respectively, +completes the set of equations for LS. + .. _methods-shannon-entropy-random-ray: ----------------------------- @@ -789,7 +953,7 @@ sources is adjusted such that: :label: fraction-source-random-ray S_i = \frac{\text{Fission source in FSR $i \times$ Volume of FSR - $i$}}{\text{Total fission source}} = \frac{Q_{i} V_{i}}{\sum_{i=1}^{i=N} + $i$}}{\text{Total fission source}} = \frac{Q_{i} V_{i}}{\sum_{i=1}^{i=N} Q_{i} V_{i}} The Shannon entropy is then computed normally as @@ -852,13 +1016,13 @@ in random ray particle transport are: areas typically have solutions that are highly effective at mitigating bias, error stemming from multigroup energy discretization is much harder to remedy. - - **Flat Source Approximation:**. In OpenMC, a "flat" (0th order) source - approximation is made, wherein the scattering and fission sources within a + - **Source Approximation:**. In OpenMC, a "flat" (0th order) source + approximation is often made, wherein the scattering and fission sources within a cell are assumed to be spatially uniform. As the source in reality is a continuous function, this leads to bias, although the bias can be reduced to acceptable levels if the flat source regions are sufficiently small. - The bias can also be mitigated by assuming a higher-order source (e.g., - linear or quadratic), although OpenMC does not yet have this capability. + The bias can also be mitigated by assuming a higher-order source such as the + linear source approximation currently implemented into OpenMC. In practical terms, this source of bias can become very large if cells are large (with dimensions beyond that of a typical particle mean free path), but the subdivision of cells can often reduce this bias to trivial levels. @@ -882,6 +1046,8 @@ in random ray particle transport are: .. _Tramm-2018: https://dspace.mit.edu/handle/1721.1/119038 .. _Tramm-2020: https://doi.org/10.1051/EPJCONF/202124703021 .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618 +.. _Ferrer-2016: https://doi.org/10.13182/NSE15-6 +.. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030 .. only:: html diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 6df132859c7..117d5e23fb5 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -447,6 +447,34 @@ in the `OpenMC Jupyter notebook collection separate materials can be defined each with a separate multigroup dataset corresponding to a given temperature. +-------------- +Linear Sources +-------------- + +Linear Sources (LS), are supported with the eigenvalue and fixed source random +ray solvers. General 3D LS can be toggled by setting the ``source_shape`` field +in the :attr:`openmc.Settings.random_ray` dictionary to ``'linear'`` as:: + + settings.random_ray['source_shape'] = 'linear' + +LS enables the use of coarser mesh discretizations and lower ray populations, +offsetting the increased computation per ray. + +While OpenMC has no specific mode for 2D simulations, such simulations can be +performed implicitly by leaving one of the dimensions of the geometry unbounded +or by imposing reflective boundary conditions with no variation in between them +in that dimension. When 3D linear sources are used in a 2D random ray +simulation, the extremely long (or potentially infinite) spatial dimension along +one of the axes can cause the linear source to become noisy, leading to +potentially large increases in variance. To mitigate this, the user can force +the z-terms of the linear source to zero by setting the ``source_shape`` field +as:: + + settings.random_ray['source_shape'] = 'linear_xy' + +which will greatly improve the quality of the linear source term in 2D +simulations. + --------------------------------- Fixed Source and Eigenvalue Modes --------------------------------- diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 1c683e5f50c..e502506c91e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -342,6 +342,8 @@ enum class RunMode { enum class SolverType { MONTE_CARLO, RANDOM_RAY }; +enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY }; + //============================================================================== // Geometry Constants diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 51938404037..33c5661dcd4 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -89,25 +89,28 @@ struct TallyTask { class FlatSourceDomain { public: //---------------------------------------------------------------------------- - // Constructors + // Constructors and Destructors FlatSourceDomain(); + virtual ~FlatSourceDomain() = default; //---------------------------------------------------------------------------- // Methods - void update_neutron_source(double k_eff); + virtual void update_neutron_source(double k_eff); double compute_k_eff(double k_eff_old) const; - void normalize_scalar_flux_and_volumes( + virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); - int64_t add_source_to_scalar_flux(); - void batch_reset(); + virtual int64_t add_source_to_scalar_flux(); + virtual void batch_reset(); void convert_source_regions_to_tallies(); void reset_tally_volumes(); void random_ray_tally(); - void accumulate_iteration_flux(); + virtual void accumulate_iteration_flux(); void output_to_vtk() const; - void all_reduce_replicated_source_regions(); + virtual void all_reduce_replicated_source_regions(); void convert_external_sources(); void count_external_source_regions(); + virtual void flux_swap(); + virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; //---------------------------------------------------------------------------- @@ -131,6 +134,7 @@ class FlatSourceDomain { vector lock_; vector was_hit_; vector volume_; + vector volume_t_; vector position_recorded_; vector position_; @@ -141,7 +145,7 @@ class FlatSourceDomain { vector source_; vector external_source_; -private: +protected: //---------------------------------------------------------------------------- // Methods void apply_external_source_to_source_region( @@ -174,7 +178,6 @@ class FlatSourceDomain { // 1D arrays representing values for all source regions vector material_; - vector volume_t_; // 2D arrays stored in 1D representing values for all source regions x energy // groups diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h new file mode 100644 index 00000000000..5010ffddd6f --- /dev/null +++ b/include/openmc/random_ray/linear_source_domain.h @@ -0,0 +1,61 @@ +#ifndef OPENMC_RANDOM_RAY_LINEAR_SOURCE_DOMAIN_H +#define OPENMC_RANDOM_RAY_LINEAR_SOURCE_DOMAIN_H + +#include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/moment_matrix.h" + +#include "openmc/openmp_interface.h" +#include "openmc/position.h" +#include "openmc/source.h" + +namespace openmc { + +/* + * The LinearSourceDomain class encompasses data and methods for storing + * scalar flux and source region for all linear source regions in a + * random ray simulation domain. + */ + +class LinearSourceDomain : public FlatSourceDomain { +public: + //---------------------------------------------------------------------------- + // Constructors + LinearSourceDomain(); + + //---------------------------------------------------------------------------- + // Methods + void update_neutron_source(double k_eff) override; + double compute_k_eff(double k_eff_old) const; + void normalize_scalar_flux_and_volumes( + double total_active_distance_per_iteration) override; + int64_t add_source_to_scalar_flux() override; + void batch_reset() override; + void convert_source_regions_to_tallies(); + void reset_tally_volumes(); + void random_ray_tally(); + void accumulate_iteration_flux() override; + void output_to_vtk() const; + void all_reduce_replicated_source_regions() override; + void convert_external_sources(); + void count_external_source_regions(); + void flux_swap() override; + double evaluate_flux_at_point(Position r, int64_t sr, int g) const override; + + //---------------------------------------------------------------------------- + // Public Data members + + vector source_gradients_; + vector flux_moments_old_; + vector flux_moments_new_; + vector flux_moments_t_; + vector centroid_; + vector centroid_iteration_; + vector centroid_t_; + vector mom_matrix_; + vector mom_matrix_t_; + +}; // class LinearSourceDomain + +} // namespace openmc + +#endif // OPENMC_RANDOM_RAY_LINEAR_SOURCE_DOMAIN_H diff --git a/include/openmc/random_ray/moment_matrix.h b/include/openmc/random_ray/moment_matrix.h new file mode 100644 index 00000000000..c95bb2c1286 --- /dev/null +++ b/include/openmc/random_ray/moment_matrix.h @@ -0,0 +1,90 @@ +#ifndef OPENMC_MOMENT_MATRIX_H +#define OPENMC_MOMENT_MATRIX_H + +#include + +#include "openmc/position.h" + +namespace openmc { + +// The MomentArray class is a 3-element array representing the x, y, and z +// moments. It is defined as an alias for the Position class to allow for +// dot products and other operations with Position objects. +// TODO: This class could in theory have 32-bit instead of 64-bit FP values. +using MomentArray = Position; + +// The MomentMatrix class is a sparse representation a 3x3 symmetric +// matrix, with elements labeled as follows: +// +// | a b c | +// | b d e | +// | c e f | +// +// This class uses FP64 values as objects that are accumulated to over many +// iterations. +class MomentMatrix { +public: + //---------------------------------------------------------------------------- + // Public data members + double a; + double b; + double c; + double d; + double e; + double f; + + //---------------------------------------------------------------------------- + // Constructors + MomentMatrix() = default; + MomentMatrix(double a, double b, double c, double d, double e, double f) + : a {a}, b {b}, c {c}, d {d}, e {e}, f {f} + {} + + //---------------------------------------------------------------------------- + // Methods + MomentMatrix inverse() const; + double determinant() const; + void compute_spatial_moments_matrix( + const Position& r, const Direction& u, const double& distance); + + inline void set_to_zero() { a = b = c = d = e = f = 0; } + + inline MomentMatrix& operator*=(double x) + { + a *= x; + b *= x; + c *= x; + d *= x; + e *= x; + f *= x; + return *this; + } + + inline MomentMatrix operator*(double x) const + { + MomentMatrix m_copy = *this; + m_copy *= x; + return m_copy; + } + + inline MomentMatrix& operator+=(const MomentMatrix& rhs) + { + a += rhs.a; + b += rhs.b; + c += rhs.c; + d += rhs.d; + e += rhs.e; + f += rhs.f; + return *this; + } + + MomentArray operator*(const MomentArray& rhs) const + { + return {a * rhs.x + b * rhs.y + c * rhs.z, + b * rhs.x + d * rhs.y + e * rhs.z, c * rhs.x + e * rhs.y + f * rhs.z}; + } +}; + +} // namespace openmc + +#endif // OPENMC_MOMENT_MATRIX_H diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index 5ee64574ec1..913a9af4a75 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -4,6 +4,7 @@ #include "openmc/memory.h" #include "openmc/particle.h" #include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/moment_matrix.h" #include "openmc/source.h" namespace openmc { @@ -25,14 +26,18 @@ class RandomRay : public Particle { // 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_linear_source(double distance, bool is_active); + void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain); uint64_t transport_history_based_single_ray(); //---------------------------------------------------------------------------- // Static data members - static double distance_inactive_; // Inactive (dead zone) ray length - static double distance_active_; // Active ray length - static unique_ptr ray_source_; // Starting source for ray sampling + static double distance_inactive_; // Inactive (dead zone) ray length + static double distance_active_; // Active ray length + static unique_ptr ray_source_; // Starting source for ray sampling + static RandomRaySourceShape source_shape_; // Flag for linear source //---------------------------------------------------------------------------- // Public data members @@ -42,6 +47,8 @@ class RandomRay : public Particle { //---------------------------------------------------------------------------- // Private data members vector delta_psi_; + vector delta_moments_; + int negroups_; FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source // data needed for ray transport diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b439574bfca..c1d47821d7a 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -2,6 +2,7 @@ #define OPENMC_RANDOM_RAY_SIMULATION_H #include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/linear_source_domain.h" namespace openmc { @@ -31,7 +32,7 @@ class RandomRaySimulation { // Data members private: // Contains all flat source region data - FlatSourceDomain domain_; + unique_ptr domain_; // Random ray eigenvalue double k_eff_ {1.0}; diff --git a/openmc/settings.py b/openmc/settings.py index 9731f54ecb1..8076625818a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -157,6 +157,9 @@ class Settings: :ray_source: Starting ray distribution (must be uniform in space and angle) as specified by a :class:`openmc.SourceBase` object. + :source_shape: + Assumed shape of the source distribution within each source + region. Options are 'flat' (default), 'linear', or 'linear_xy'. :volume_normalized_flux_tallies: Whether to normalize flux tallies by volume (bool). The default is 'False'. When enabled, flux tallies will be reported in units of @@ -1090,6 +1093,9 @@ def random_ray(self, random_ray: dict): random_ray[key], 0.0, True) elif key == 'ray_source': cv.check_type('random ray source', random_ray[key], SourceBase) + elif key == 'source_shape': + cv.check_value('source shape', random_ray[key], + ('flat', 'linear', 'linear_xy')) elif key == 'volume_normalized_flux_tallies': cv.check_type('volume normalized flux tallies', random_ray[key], bool) else: @@ -1885,6 +1891,8 @@ def _random_ray_from_xml_element(self, root): elif child.tag == 'source': source = SourceBase.from_xml_element(child) self.random_ray['ray_source'] = source + elif child.tag == 'source_shape': + self.random_ray['source_shape'] = child.text elif child.tag == 'volume_normalized_flux_tallies': self.random_ray['volume_normalized_flux_tallies'] = ( child.text in ('true', '1') diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 29b72c77bed..8b6cda93f2a 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -751,6 +751,13 @@ void FlatSourceDomain::all_reduce_replicated_source_regions() #endif } +double FlatSourceDomain::evaluate_flux_at_point( + Position r, int64_t sr, int g) const +{ + return scalar_flux_final_[sr * negroups_ + g] / + (settings::n_batches - settings::n_inactive); +} + // Outputs all basic material, FSR ID, multigroup flux, and // fission source data to .vtk file that can be directly // loaded and displayed by Paraview. Note that .vtk binary @@ -811,6 +818,7 @@ void FlatSourceDomain::output_to_vtk() const // Relate voxel spatial locations to random ray source regions vector voxel_indices(Nx * Ny * Nz); + vector voxel_positions(Nx * Ny * Nz); #pragma omp parallel for collapse(3) for (int z = 0; z < Nz; z++) { @@ -827,6 +835,7 @@ void FlatSourceDomain::output_to_vtk() const int64_t source_region_idx = source_region_offsets_[i_cell] + p.cell_instance(); voxel_indices[z * Ny * Nx + y * Nx + x] = source_region_idx; + voxel_positions[z * Ny * Nx + y * Nx + x] = sample; } } } @@ -851,11 +860,10 @@ void FlatSourceDomain::output_to_vtk() const for (int g = 0; g < negroups_; g++) { std::fprintf(plot, "SCALARS flux_group_%d float\n", g); std::fprintf(plot, "LOOKUP_TABLE default\n"); - for (int fsr : voxel_indices) { + for (int i = 0; i < Nx * Ny * Nz; i++) { + int64_t fsr = voxel_indices[i]; int64_t source_element = fsr * negroups_ + g; - float flux = - scalar_flux_final_[source_element] * source_normalization_factor; - flux /= (settings::n_batches - settings::n_inactive); + float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g); flux = convert_to_big_endian(flux); std::fwrite(&flux, sizeof(float), 1, plot); } @@ -882,14 +890,14 @@ void FlatSourceDomain::output_to_vtk() const // Plot fission source std::fprintf(plot, "SCALARS total_fission_source float\n"); std::fprintf(plot, "LOOKUP_TABLE default\n"); - for (int fsr : voxel_indices) { + for (int i = 0; i < Nx * Ny * Nz; i++) { + int64_t fsr = voxel_indices[i]; + float total_fission = 0.0; int mat = material_[fsr]; for (int g = 0; g < negroups_; g++) { int64_t source_element = fsr * negroups_ + g; - float flux = - scalar_flux_final_[source_element] * source_normalization_factor; - flux /= (settings::n_batches - settings::n_inactive); + float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g); float Sigma_f = data::mg.macro_xs_[mat].get_xs( MgxsType::FISSION, g, nullptr, nullptr, nullptr, 0, 0); total_fission += Sigma_f * flux; @@ -1027,5 +1035,9 @@ void FlatSourceDomain::convert_external_sources() } } } +void FlatSourceDomain::flux_swap() +{ + scalar_flux_old_.swap(scalar_flux_new_); +} } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp new file mode 100644 index 00000000000..1603ec24a4e --- /dev/null +++ b/src/random_ray/linear_source_domain.cpp @@ -0,0 +1,269 @@ +#include "openmc/random_ray/linear_source_domain.h" + +#include "openmc/cell.h" +#include "openmc/geometry.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" +#include "openmc/output.h" +#include "openmc/plot.h" +#include "openmc/random_ray/random_ray.h" +#include "openmc/simulation.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" +#include "openmc/timer.h" + +namespace openmc { + +//============================================================================== +// LinearSourceDomain implementation +//============================================================================== + +LinearSourceDomain::LinearSourceDomain() : FlatSourceDomain() +{ + // First order spatial moment of the scalar flux + flux_moments_old_.assign(n_source_elements_, {0.0, 0.0, 0.0}); + flux_moments_new_.assign(n_source_elements_, {0.0, 0.0, 0.0}); + flux_moments_t_.assign(n_source_elements_, {0.0, 0.0, 0.0}); + // Source gradients given by M inverse multiplied by source moments + source_gradients_.assign(n_source_elements_, {0.0, 0.0, 0.0}); + + centroid_.assign(n_source_regions_, {0.0, 0.0, 0.0}); + centroid_iteration_.assign(n_source_regions_, {0.0, 0.0, 0.0}); + centroid_t_.assign(n_source_regions_, {0.0, 0.0, 0.0}); + mom_matrix_.assign(n_source_regions_, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); + mom_matrix_t_.assign(n_source_regions_, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); +} + +void LinearSourceDomain::batch_reset() +{ + FlatSourceDomain::batch_reset(); +#pragma omp parallel for + for (int64_t se = 0; se < n_source_elements_; se++) { + flux_moments_new_[se] = {0.0, 0.0, 0.0}; + } +#pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + centroid_iteration_[sr] = {0.0, 0.0, 0.0}; + mom_matrix_[sr] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + } +} + +void LinearSourceDomain::update_neutron_source(double k_eff) +{ + simulation::time_update_src.start(); + + double inverse_k_eff = 1.0 / k_eff; + + // Temperature and angle indices, if using multiple temperature + // data sets and/or anisotropic data sets. + // TODO: Currently assumes we are only using single temp/single + // angle data. + const int t = 0; + const int a = 0; + +#pragma omp parallel for + for (int sr = 0; sr < n_source_regions_; sr++) { + + int material = material_[sr]; + MomentMatrix invM = mom_matrix_[sr].inverse(); + + for (int e_out = 0; e_out < negroups_; e_out++) { + float sigma_t = data::mg.macro_xs_[material].get_xs( + MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a); + + float scatter_flat = 0.0f; + float fission_flat = 0.0f; + MomentArray scatter_linear = {0.0, 0.0, 0.0}; + MomentArray fission_linear = {0.0, 0.0, 0.0}; + + for (int e_in = 0; e_in < negroups_; e_in++) { + // Handles for the flat and linear components of the flux + float flux_flat = scalar_flux_old_[sr * negroups_ + e_in]; + MomentArray flux_linear = flux_moments_old_[sr * negroups_ + e_in]; + + // Handles for cross sections + float sigma_s = data::mg.macro_xs_[material].get_xs( + MgxsType::NU_SCATTER, e_in, &e_out, nullptr, nullptr, t, a); + float nu_sigma_f = data::mg.macro_xs_[material].get_xs( + MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a); + float chi = data::mg.macro_xs_[material].get_xs( + MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a); + + // Compute source terms for flat and linear components of the flux + scatter_flat += sigma_s * flux_flat; + fission_flat += nu_sigma_f * flux_flat * chi; + scatter_linear += sigma_s * flux_linear; + fission_linear += nu_sigma_f * flux_linear * chi; + } + + // Compute the flat source term + source_[sr * negroups_ + e_out] = + (scatter_flat + fission_flat * inverse_k_eff) / sigma_t; + + // Compute the linear source terms + if (simulation::current_batch > 2) { + source_gradients_[sr * negroups_ + e_out] = + invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); + } + } + } + + if (settings::run_mode == RunMode::FIXED_SOURCE) { +// Add external source to flat source term if in fixed source mode +#pragma omp parallel for + for (int se = 0; se < n_source_elements_; se++) { + source_[se] += external_source_[se]; + } + } + + simulation::time_update_src.stop(); +} + +void LinearSourceDomain::normalize_scalar_flux_and_volumes( + double total_active_distance_per_iteration) +{ + float normalization_factor = 1.0 / total_active_distance_per_iteration; + double volume_normalization_factor = + 1.0 / (total_active_distance_per_iteration * simulation::current_batch); + +// Normalize flux to total distance travelled by all rays this iteration +#pragma omp parallel for + for (int64_t e = 0; e < scalar_flux_new_.size(); e++) { + scalar_flux_new_[e] *= normalization_factor; + flux_moments_new_[e] *= 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++) { + centroid_t_[sr] += centroid_iteration_[sr]; + mom_matrix_t_[sr] += mom_matrix_[sr]; + volume_t_[sr] += volume_[sr]; + volume_[sr] = volume_t_[sr] * volume_normalization_factor; + if (volume_t_[sr] > 0.0) { + double inv_volume = 1.0 / volume_t_[sr]; + centroid_[sr] = centroid_t_[sr]; + centroid_[sr] *= inv_volume; + mom_matrix_[sr] = mom_matrix_t_[sr]; + mom_matrix_[sr] *= inv_volume; + } + } +} + +int64_t LinearSourceDomain::add_source_to_scalar_flux() +{ + int64_t n_hits = 0; + + // Temperature and angle indices, if using multiple temperature + // data sets and/or anisotropic data sets. + // TODO: Currently assumes we are only using single temp/single + // angle data. + const int t = 0; + const int a = 0; + +#pragma omp parallel for reduction(+ : n_hits) + for (int sr = 0; sr < n_source_regions_; sr++) { + + double volume = volume_[sr]; + int material = material_[sr]; + + // Check if this cell was hit this iteration + int was_cell_hit = was_hit_[sr]; + if (was_cell_hit) { + n_hits++; + } + + for (int g = 0; g < negroups_; g++) { + int64_t idx = (sr * negroups_) + g; + // There are three scenarios we need to consider: + if (was_cell_hit) { + // 1. If the FSR was hit this iteration, then the new flux is equal to + // the flat source from the previous iteration plus the contributions + // from rays passing through the source region (computed during the + // transport sweep) + scalar_flux_new_[idx] /= volume; + scalar_flux_new_[idx] += source_[idx]; + flux_moments_new_[idx] *= (1.0 / volume); + } else if (volume > 0.0) { + // 2. If the FSR was not hit this iteration, but has been hit some + // previous iteration, then we simply set the new scalar flux to be + // equal to the contribution from the flat source alone. + scalar_flux_new_[idx] = source_[idx]; + } else { + // 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. + scalar_flux_new_[idx] = 0.0f; + flux_moments_new_[idx] *= 0.0; + } + } + } + + return n_hits; +} + +void LinearSourceDomain::flux_swap() +{ + FlatSourceDomain::flux_swap(); + flux_moments_old_.swap(flux_moments_new_); +} + +void LinearSourceDomain::accumulate_iteration_flux() +{ + // Accumulate scalar flux + FlatSourceDomain::accumulate_iteration_flux(); + + // Accumulate scalar flux moments +#pragma omp parallel for + for (int64_t se = 0; se < n_source_elements_; se++) { + flux_moments_t_[se] += flux_moments_new_[se]; + } +} + +void LinearSourceDomain::all_reduce_replicated_source_regions() +{ +#ifdef OPENMC_MPI + FlatSourceDomain::all_reduce_replicated_source_regions(); + simulation::time_bank_sendrecv.start(); + + // 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(flux_moments_new_.data()), + n_source_elements_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce(MPI_IN_PLACE, static_cast(mom_matrix_.data()), + n_source_regions_ * 6, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce(MPI_IN_PLACE, static_cast(centroid_iteration_.data()), + n_source_regions_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + + simulation::time_bank_sendrecv.stop(); +#endif +} + +double LinearSourceDomain::evaluate_flux_at_point( + Position r, int64_t sr, int g) const +{ + float phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g); + + Position local_r = r - centroid_[sr]; + MomentArray phi_linear = flux_moments_t_[sr * negroups_ + g]; + phi_linear *= 1.0 / (settings::n_batches - settings::n_inactive); + + MomentMatrix invM = mom_matrix_[sr].inverse(); + MomentArray phi_solved = invM * phi_linear; + + return phi_flat + phi_solved.dot(local_r); +} + +} // namespace openmc diff --git a/src/random_ray/moment_matrix.cpp b/src/random_ray/moment_matrix.cpp new file mode 100644 index 00000000000..0324a14943b --- /dev/null +++ b/src/random_ray/moment_matrix.cpp @@ -0,0 +1,84 @@ +#include "openmc/random_ray/moment_matrix.h" +#include "openmc/error.h" + +#include + +namespace openmc { + +//============================================================================== +// UpperTriangular implementation +//============================================================================== + +// Inverts a 3x3 smmetric matrix labeled as: +// +// | a b c | +// | b d e | +// | c e f | +// +// We first check the determinant to ensure it is non-zero before proceeding +// with the inversion. If the determinant is zero, we return a matrix of zeros. +// Inversion is calculated by computing the adjoint matrix first, and then the +// inverse can be computed as: A^-1 = 1/det(A) * adj(A) +MomentMatrix MomentMatrix::inverse() const +{ + MomentMatrix inv; + + // Check if the determinant is zero + double det = determinant(); + if (det < std::abs(1.0e-10)) { + // Set the inverse to zero. In effect, this will + // result in all the linear terms of the source becoming + // zero, leaving just the flat source. + inv.set_to_zero(); + return inv; + } + + // Compute the adjoint matrix + inv.a = d * f - e * e; + inv.b = c * e - b * f; + inv.c = b * e - c * d; + inv.d = a * f - c * c; + inv.e = b * c - a * e; + inv.f = a * d - b * b; + + // A^-1 = 1/det(A) * adj(A) + inv *= 1.0 / det; + + return inv; +} + +// Computes the determinant of a 3x3 symmetric +// matrix, with elements labeled as follows: +// +// | a b c | +// | b d e | +// | c e f | +double MomentMatrix::determinant() const +{ + return a * (d * f - e * e) - b * (b * f - c * e) + c * (b * e - c * d); +} + +// Compute a 3x3 spatial moment matrix based on a single ray crossing. +// The matrix is symmetric, and is defined as: +// +// | a b c | +// | b d e | +// | c e f | +// +// The estimate of the obect's spatial moments matrix is computed based on the +// midpoint of the ray's crossing, the direction of the ray, and the distance +// the ray traveled through the 3D object. +void MomentMatrix::compute_spatial_moments_matrix( + const Position& r, const Direction& u, const double& distance) +{ + constexpr double one_over_twelve = 1.0 / 12.0; + const double distance2_12 = distance * distance * one_over_twelve; + a = r[0] * r[0] + u[0] * u[0] * distance2_12; + b = r[0] * r[1] + u[0] * u[1] * distance2_12; + c = r[0] * r[2] + u[0] * u[2] * distance2_12; + d = r[1] * r[1] + u[1] * u[1] * distance2_12; + e = r[1] * r[2] + u[1] * u[2] * distance2_12; + f = r[2] * r[2] + u[2] * u[2] * distance2_12; +} + +} // namespace openmc \ No newline at end of file diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 63d728cce8b..a5bf6ec1060 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -1,9 +1,11 @@ #include "openmc/random_ray/random_ray.h" +#include "openmc/constants.h" #include "openmc/geometry.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/linear_source_domain.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" @@ -60,6 +62,118 @@ float cjosey_exponential(float tau) return num / den; } +// The below two functions (exponentialG and exponentialG2) were developed +// by Colin Josey. The implementation of these functions is closely based +// on the OpenMOC versions of these functions. The OpenMOC license is given +// below: + +// Copyright (C) 2012-2023 Massachusetts Institute of Technology and OpenMOC +// contributors +// +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// +// The above copyright notice and this permission notice shall be included in +// all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +// SOFTWARE. + +// Computes y = 1/x-(1-exp(-x))/x**2 using a 5/6th order rational +// approximation. It is accurate to 2e-7 over [0, 1e5]. Developed by Colin +// Josey using Remez's algorithm, with original implementation in OpenMOC at: +// https://github.com/mit-crpg/OpenMOC/blob/develop/src/exponentials.h +float exponentialG(float tau) +{ + // Numerator coefficients in rational approximation for 1/x - (1 - exp(-x)) / + // x^2 + constexpr float d0n = 0.5f; + constexpr float d1n = 0.176558112351595f; + constexpr float d2n = 0.04041584305811143f; + constexpr float d3n = 0.006178333902037397f; + constexpr float d4n = 0.0006429894635552992f; + constexpr float d5n = 0.00006064409107557148f; + + // Denominator coefficients in rational approximation for 1/x - (1 - exp(-x)) + // / x^2 + constexpr float d0d = 1.0f; + constexpr float d1d = 0.6864462055546078f; + constexpr float d2d = 0.2263358514260129f; + constexpr float d3d = 0.04721469893686252f; + constexpr float d4d = 0.006883236664917246f; + constexpr float d5d = 0.0007036272419147752f; + constexpr float d6d = 0.00006064409107557148f; + + float x = tau; + + float num = d5n; + num = num * x + d4n; + num = num * x + d3n; + num = num * x + d2n; + num = num * x + d1n; + num = num * x + d0n; + + float den = d6d; + den = den * x + d5d; + den = den * x + d4d; + den = den * x + d3d; + den = den * x + d2d; + den = den * x + d1d; + den = den * x + d0d; + + return num / den; +} + +// Computes G2 : y = 2/3 - (1 + 2/x) * (1/x + 0.5 - (1 + 1/x) * (1-exp(-x)) / +// x) using a 5/5th order rational approximation. It is accurate to 1e-6 over +// [0, 1e6]. Developed by Colin Josey using Remez's algorithm, with original +// implementation in OpenMOC at: +// https://github.com/mit-crpg/OpenMOC/blob/develop/src/exponentials.h +float exponentialG2(float tau) +{ + + // Coefficients for numerator in rational approximation + constexpr float g1n = -0.08335775885589858f; + constexpr float g2n = -0.003603942303847604f; + constexpr float g3n = 0.0037673183263550827f; + constexpr float g4n = 0.00001124183494990467f; + constexpr float g5n = 0.00016837426505799449f; + + // Coefficients for denominator in rational approximation + constexpr float g1d = 0.7454048371823628f; + constexpr float g2d = 0.23794300531408347f; + constexpr float g3d = 0.05367250964303789f; + constexpr float g4d = 0.006125197988351906f; + constexpr float g5d = 0.0010102514456857377f; + + float x = tau; + + float num = g5n; + num = num * x + g4n; + num = num * x + g3n; + num = num * x + g2n; + num = num * x + g1n; + num = num * x; + + float den = g5d; + den = den * x + g4d; + den = den * x + g3d; + den = den * x + g2d; + den = den * x + g1d; + den = den * x + 1.0f; + + return num / den; +} + //============================================================================== // RandomRay implementation //============================================================================== @@ -68,12 +182,18 @@ float cjosey_exponential(float tau) double RandomRay::distance_inactive_; double RandomRay::distance_active_; unique_ptr RandomRay::ray_source_; +RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT}; RandomRay::RandomRay() : angular_flux_(data::mg.num_energy_groups_), delta_psi_(data::mg.num_energy_groups_), negroups_(data::mg.num_energy_groups_) -{} +{ + if (source_shape_ == RandomRaySourceShape::LINEAR || + source_shape_ == RandomRaySourceShape::LINEAR_XY) { + delta_moments_.resize(negroups_); + } +} RandomRay::RandomRay(uint64_t ray_id, FlatSourceDomain* domain) : RandomRay() { @@ -152,6 +272,21 @@ void RandomRay::event_advance_ray() } } +void RandomRay::attenuate_flux(double distance, bool is_active) +{ + switch (source_shape_) { + case RandomRaySourceShape::FLAT: + attenuate_flux_flat_source(distance, is_active); + break; + case RandomRaySourceShape::LINEAR: + case RandomRaySourceShape::LINEAR_XY: + attenuate_flux_linear_source(distance, is_active); + break; + default: + fatal_error("Unknown source shape for random ray transport."); + } +} + // This function forms the inner loop of the random ray transport process. // It is responsible for several tasks. Based on the incoming angular flux // of the ray and the source term in the region, the outgoing angular flux @@ -165,7 +300,7 @@ void RandomRay::event_advance_ray() // 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(double distance, bool is_active) +void RandomRay::attenuate_flux_flat_source(double distance, bool is_active) { // The number of geometric intersections is counted for reporting purposes n_event()++; @@ -236,6 +371,158 @@ void RandomRay::attenuate_flux(double distance, bool is_active) } } +void RandomRay::attenuate_flux_linear_source(double distance, bool is_active) +{ + // Cast domain to LinearSourceDomain + LinearSourceDomain* domain = dynamic_cast(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 source_region = + domain_->source_region_offsets_[i_cell] + cell_instance(); + + // The source element is the energy-specific region index + int64_t source_element = source_region * negroups_; + int material = this->material(); + + // Temperature and angle indices, if using multiple temperature + // data sets and/or anisotropic data sets. + // TODO: Currently assumes we are only using single temp/single + // angle data. + const int t = 0; + const int a = 0; + + Position& centroid = domain->centroid_[source_region]; + 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 + Position rm_local; + Position r0_local; + + // In the first few iterations of the simulation, the source region + // may not yet have had any ray crossings, in which case there will + // 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->volume_t_[source_region]) { + rm_local = midpoint - centroid; + r0_local = r() - centroid; + } else { + rm_local = {0.0, 0.0, 0.0}; + r0_local = -u() * 0.5 * distance; + } + double distance_2 = distance * distance; + + // Linear Source MOC incoming flux attenuation + source + // contribution/attenuation equation + for (int g = 0; g < negroups_; g++) { + + // Compute tau, the optical thickness of the ray segment + float sigma_t = data::mg.macro_xs_[material].get_xs( + MgxsType::TOTAL, g, NULL, NULL, NULL, t, a); + float tau = sigma_t * distance; + + // If tau is very small, set it to zero to avoid numerical issues. + // The following computations will still work with tau = 0. + if (tau < 1.0e-8f) { + tau = 0.0f; + } + + // Compute linear source terms, spatial and directional (dir), + // calculated from the source gradients dot product with local centroid + // and direction, respectively. + float spatial_source = + domain_->source_[source_element + g] + + rm_local.dot(domain->source_gradients_[source_element + g]); + float dir_source = u().dot(domain->source_gradients_[source_element + g]); + + float gn = exponentialG(tau); + float f1 = 1.0f - tau * gn; + float f2 = (2.0f * gn - f1) * distance_2; + float new_delta_psi = (angular_flux_[g] - spatial_source) * f1 * distance - + 0.5 * dir_source * f2; + + float h1 = f1 - gn; + float g1 = 0.5f - h1; + float g2 = exponentialG2(tau); + g1 = g1 * spatial_source; + g2 = g2 * dir_source * distance * 0.5f; + h1 = h1 * angular_flux_[g]; + h1 = (g1 + g2 + h1) * distance_2; + spatial_source = spatial_source * distance + new_delta_psi; + + // Store contributions for this group into arrays, so that they can + // be accumulated into the source region's estimates inside of the locked + // region. + delta_psi_[g] = new_delta_psi; + delta_moments_[g] = r0_local * spatial_source + u() * h1; + + // Update the angular flux for this group + angular_flux_[g] -= new_delta_psi * sigma_t; + + // If 2D mode is enabled, the z-component of the flux moments is forced + // to zero + if (source_shape_ == RandomRaySourceShape::LINEAR_XY) { + delta_moments_[g].z = 0.0; + } + } + + // 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_->lock_[source_region].lock(); + + // Accumulate deltas into the new estimate of source region flux for this + // iteration + for (int g = 0; g < negroups_; g++) { + domain_->scalar_flux_new_[source_element + g] += delta_psi_[g]; + domain->flux_moments_new_[source_element + 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_->volume_[source_region] += distance; + domain->centroid_iteration_[source_region] += midpoint * distance; + moment_matrix_estimate *= distance; + domain->mom_matrix_[source_region] += moment_matrix_estimate; + + // If the source region hasn't been hit yet this iteration, + // indicate that it now has + if (domain_->was_hit_[source_region] == 0) { + domain_->was_hit_[source_region] = 1; + } + + // Tally valid position inside the source region (e.g., midpoint of + // the ray) if not done already + if (!domain_->position_recorded_[source_region]) { + domain_->position_[source_region] = midpoint; + domain_->position_recorded_[source_region] = 1; + } + + // Release lock + domain_->lock_[source_region].unlock(); + } +} + void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) { domain_ = domain; diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 17866c30208..4bc77645bcd 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -6,6 +6,7 @@ #include "openmc/mgxs_interface.h" #include "openmc/output.h" #include "openmc/plot.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/random_ray/random_ray.h" #include "openmc/simulation.h" #include "openmc/source.h" @@ -241,14 +242,26 @@ RandomRaySimulation::RandomRaySimulation() // Random ray mode does not have an inner loop over generations within a // batch, so set the current gen to 1 simulation::current_gen = 1; + + switch (RandomRay::source_shape_) { + case RandomRaySourceShape::FLAT: + domain_ = make_unique(); + break; + case RandomRaySourceShape::LINEAR: + case RandomRaySourceShape::LINEAR_XY: + domain_ = make_unique(); + break; + default: + fatal_error("Unknown random ray source shape"); + } } void RandomRaySimulation::simulate() { 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_->convert_external_sources(); + domain_->count_external_source_regions(); } // Random ray power iteration loop @@ -262,11 +275,11 @@ void RandomRaySimulation::simulate() simulation::total_weight = 1.0; // Update source term (scattering + fission) - domain_.update_neutron_source(k_eff_); + domain_->update_neutron_source(k_eff_); // Reset scalar fluxes, iteration volume tallies, and region hit flags to // zero - domain_.batch_reset(); + domain_->batch_reset(); // Start timer for transport simulation::time_transport.start(); @@ -275,7 +288,7 @@ void RandomRaySimulation::simulate() #pragma omp parallel for schedule(dynamic) \ reduction(+ : total_geometric_intersections_) for (int i = 0; i < simulation::work_per_rank; i++) { - RandomRay ray(i, &domain_); + RandomRay ray(i, domain_.get()); total_geometric_intersections_ += ray.transport_history_based_single_ray(); } @@ -283,18 +296,18 @@ void RandomRaySimulation::simulate() simulation::time_transport.stop(); // If using multiple MPI ranks, perform all reduce on all transport results - domain_.all_reduce_replicated_source_regions(); + domain_->all_reduce_replicated_source_regions(); // Normalize scalar flux and update volumes - domain_.normalize_scalar_flux_and_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(); + 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_); + 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_; @@ -304,19 +317,19 @@ void RandomRaySimulation::simulate() if (simulation::current_batch > settings::n_inactive && mpi::master) { // Generate mapping between source regions and tallies - if (!domain_.mapped_all_tallies_) { - domain_.convert_source_regions_to_tallies(); + if (!domain_->mapped_all_tallies_) { + domain_->convert_source_regions_to_tallies(); } // Use above mapping to contribute FSR flux data to appropriate tallies - domain_.random_ray_tally(); + domain_->random_ray_tally(); // Add this iteration's scalar flux estimate to final accumulated estimate - domain_.accumulate_iteration_flux(); + domain_->accumulate_iteration_flux(); } // Set phi_old = phi_new - domain_.scalar_flux_old_.swap(domain_.scalar_flux_new_); + domain_->flux_swap(); // Check for any obvious insabilities/nans/infs instability_check(n_hits, k_eff_, avg_miss_rate_); @@ -347,9 +360,9 @@ void RandomRaySimulation::output_simulation_results() const 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(); + domain_->output_to_vtk(); } } } @@ -359,25 +372,28 @@ 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(domain_.n_source_regions_)) * + double percent_missed = ((domain_->n_source_regions_ - n_hits) / + static_cast(domain_->n_source_regions_)) * 100.0; avg_miss_rate += percent_missed; - if (percent_missed > 10.0) { - warning(fmt::format( - "Very high FSR miss rate detected ({:.3f}%). Instability may occur. " - "Increase ray density by adding more rays and/or active distance.", - percent_missed)); - } else if (percent_missed > 0.01) { - warning(fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing " - "ray density by adding more rays and/or active " - "distance may improve simulation efficiency.", - percent_missed)); - } + if (mpi::master) { + if (percent_missed > 10.0) { + warning(fmt::format( + "Very high FSR miss rate detected ({:.3f}%). Instability may occur. " + "Increase ray density by adding more rays and/or active distance.", + percent_missed)); + } else if (percent_missed > 0.01) { + warning( + fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing " + "ray density by adding more rays and/or active " + "distance may improve simulation efficiency.", + percent_missed)); + } - if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) { - fatal_error("Instability detected"); + if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) { + fatal_error("Instability detected"); + } } } diff --git a/src/settings.cpp b/src/settings.cpp index b92b77df937..db956abf5ca 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -269,6 +269,19 @@ void get_run_parameters(pugi::xml_node node_base) } else { fatal_error("Specify random ray source in settings XML"); } + if (check_for_node(random_ray_node, "source_shape")) { + std::string temp_str = + get_node_value(random_ray_node, "source_shape", true, true); + if (temp_str == "flat") { + RandomRay::source_shape_ = RandomRaySourceShape::FLAT; + } else if (temp_str == "linear") { + RandomRay::source_shape_ = RandomRaySourceShape::LINEAR; + } else if (temp_str == "linear_xy") { + RandomRay::source_shape_ = RandomRaySourceShape::LINEAR_XY; + } else { + fatal_error("Unrecognized source shape: " + temp_str); + } + } if (check_for_node(random_ray_node, "volume_normalized_flux_tallies")) { FlatSourceDomain::volume_normalized_flux_tallies_ = get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies"); diff --git a/tests/regression_tests/random_ray_fixed_source_linear/__init__.py b/tests/regression_tests/random_ray_fixed_source_linear/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat new file mode 100644 index 00000000000..e2972c5f20a --- /dev/null +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear/inputs_true.dat @@ -0,0 +1,245 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +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 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + True + linear + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear/results_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear/results_true.dat new file mode 100644 index 00000000000..7fb01a0ea11 --- /dev/null +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +-5.745718E+02 +9.757661E+04 +tally 2: +3.074428E-02 +1.952251E-04 +tally 3: +1.980876E-03 +7.970502E-07 diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat new file mode 100644 index 00000000000..5b958d65cdb --- /dev/null +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/inputs_true.dat @@ -0,0 +1,245 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +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 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + True + linear_xy + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/results_true.dat b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/results_true.dat new file mode 100644 index 00000000000..22b39edcfd8 --- /dev/null +++ b/tests/regression_tests/random_ray_fixed_source_linear/linear_xy/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +-5.745810E+02 +9.758220E+04 +tally 2: +3.022777E-02 +1.884091E-04 +tally 3: +1.980651E-03 +7.968779E-07 diff --git a/tests/regression_tests/random_ray_fixed_source_linear/test.py b/tests/regression_tests/random_ray_fixed_source_linear/test.py new file mode 100644 index 00000000000..25335dea602 --- /dev/null +++ b/tests/regression_tests/random_ray_fixed_source_linear/test.py @@ -0,0 +1,27 @@ +import os + +import numpy as np +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 TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +@pytest.mark.parametrize("shape", ["linear", "linear_xy"]) +def test_random_ray_fixed_source_linear(shape): + with change_directory(shape): + openmc.reset_auto_ids() + model = random_ray_three_region_cube() + model.settings.random_ray['source_shape'] = shape + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/random_ray_linear/__init__.py b/tests/regression_tests/random_ray_linear/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_linear/linear/inputs_true.dat b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat new file mode 100644 index 00000000000..c580be3d19e --- /dev/null +++ b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat @@ -0,0 +1,110 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.126 0.126 + 10 10 + -0.63 -0.63 + +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 + + + 1.26 1.26 + 2 2 + -1.26 -1.26 + +2 2 +2 5 + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + multi-group + + 100.0 + 20.0 + + + -1.26 -1.26 -1 1.26 1.26 1 + + + True + linear + + + + + 2 2 + -1.26 -1.26 + 1.26 1.26 + + + 1 + + + 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + + + 1 2 + flux fission nu-fission + analog + + + diff --git a/tests/regression_tests/random_ray_linear/linear/results_true.dat b/tests/regression_tests/random_ray_linear/linear/results_true.dat new file mode 100644 index 00000000000..a3e24dc75d8 --- /dev/null +++ b/tests/regression_tests/random_ray_linear/linear/results_true.dat @@ -0,0 +1,171 @@ +k-combined: +8.273022E-01 1.347623E-02 +tally 1: +5.004109E+00 +5.022655E+00 +1.844047E+00 +6.833376E-01 +4.488042E+00 +4.047669E+00 +2.824818E+00 +1.599927E+00 +4.182704E-01 +3.509796E-02 +1.017987E+00 +2.078987E-01 +1.676761E+00 +5.682729E-01 +5.385624E-02 +5.857594E-04 +1.310753E-01 +3.469677E-03 +2.353602E+00 +1.127695E+00 +7.721312E-02 +1.212202E-03 +1.879213E-01 +7.180339E-03 +7.082957E+00 +1.019023E+01 +8.203580E-02 +1.366996E-03 +1.996612E-01 +8.097436E-03 +2.034293E+01 +8.321967E+01 +3.099116E-02 +1.933713E-04 +7.668546E-02 +1.183975E-03 +1.311478E+01 +3.442575E+01 +1.778200E-01 +6.347187E-03 +4.945973E-01 +4.910476E-02 +7.576546E+00 +1.148094E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.386885E+00 +2.295318E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.823720E+00 +6.769694E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.703042E+00 +1.494228E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +7.465612E+00 +1.132769E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.820337E+01 +6.657534E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.127686E+01 +2.549122E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.561053E+00 +4.169936E+00 +1.674166E+00 +5.624648E-01 +4.074585E+00 +3.331694E+00 +2.723008E+00 +1.487311E+00 +4.059013E-01 +3.307772E-02 +9.878827E-01 +1.959320E-01 +1.663524E+00 +5.584734E-01 +5.399658E-02 +5.877094E-04 +1.314169E-01 +3.481227E-03 +2.310994E+00 +1.086435E+00 +7.627297E-02 +1.181877E-03 +1.856331E-01 +7.000707E-03 +7.100201E+00 +1.023633E+01 +8.285591E-02 +1.393520E-03 +2.016572E-01 +8.254554E-03 +2.119292E+01 +9.022917E+01 +3.281214E-02 +2.163583E-04 +8.119135E-02 +1.324719E-03 +1.380975E+01 +3.815787E+01 +1.919186E-01 +7.384267E-03 +5.338118E-01 +5.712809E-02 +5.024478E+00 +5.056216E+00 +1.888826E+00 +7.148091E-01 +4.597025E+00 +4.234087E+00 +2.839930E+00 +1.616473E+00 +4.294394E-01 +3.697612E-02 +1.045170E+00 +2.190237E-01 +1.688196E+00 +5.765397E-01 +5.539650E-02 +6.200086E-04 +1.348240E-01 +3.672547E-03 +2.361633E+00 +1.136540E+00 +7.901052E-02 +1.270359E-03 +1.922958E-01 +7.524822E-03 +7.100780E+00 +1.024481E+01 +8.389610E-02 +1.429594E-03 +2.041888E-01 +8.468240E-03 +2.049187E+01 +8.438141E+01 +3.195576E-02 +2.052140E-04 +7.907229E-02 +1.256485E-03 +1.327412E+01 +3.524466E+01 +1.850084E-01 +6.847675E-03 +5.145915E-01 +5.297677E-02 diff --git a/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat new file mode 100644 index 00000000000..5c4076de70f --- /dev/null +++ b/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat @@ -0,0 +1,110 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.126 0.126 + 10 10 + -0.63 -0.63 + +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 + + + 1.26 1.26 + 2 2 + -1.26 -1.26 + +2 2 +2 5 + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + multi-group + + 100.0 + 20.0 + + + -1.26 -1.26 -1 1.26 1.26 1 + + + True + linear_xy + + + + + 2 2 + -1.26 -1.26 + 1.26 1.26 + + + 1 + + + 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + + + 1 2 + flux fission nu-fission + analog + + + diff --git a/tests/regression_tests/random_ray_linear/linear_xy/results_true.dat b/tests/regression_tests/random_ray_linear/linear_xy/results_true.dat new file mode 100644 index 00000000000..f79b3fa787f --- /dev/null +++ b/tests/regression_tests/random_ray_linear/linear_xy/results_true.dat @@ -0,0 +1,171 @@ +k-combined: +8.368882E-01 8.107070E-03 +tally 1: +5.072700E+00 +5.152684E+00 +1.876680E+00 +7.051697E-01 +4.567463E+00 +4.176989E+00 +2.858196E+00 +1.636775E+00 +4.239946E-01 +3.601884E-02 +1.031918E+00 +2.133534E-01 +1.692006E+00 +5.789664E-01 +5.442292E-02 +5.988675E-04 +1.324545E-01 +3.547321E-03 +2.371801E+00 +1.146513E+00 +7.804414E-02 +1.241074E-03 +1.899438E-01 +7.351359E-03 +7.134037E+00 +1.034554E+01 +8.269855E-02 +1.390936E-03 +2.012742E-01 +8.239245E-03 +2.043174E+01 +8.386878E+01 +3.098171E-02 +1.929186E-04 +7.666208E-02 +1.181203E-03 +1.312800E+01 +3.447362E+01 +1.763141E-01 +6.217438E-03 +4.904088E-01 +4.810096E-02 +7.608165E+00 +1.157716E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.392187E+00 +2.302667E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.821339E+00 +6.737984E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.695742E+00 +1.483104E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +7.457904E+00 +1.129343E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.824331E+01 +6.687165E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.138096E+01 +2.591161E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.582534E+00 +4.206969E+00 +1.707129E+00 +5.837143E-01 +4.154811E+00 +3.457563E+00 +2.727459E+00 +1.491183E+00 +4.100516E-01 +3.371280E-02 +9.979836E-01 +1.996938E-01 +1.660765E+00 +5.570066E-01 +5.425719E-02 +5.941503E-04 +1.320511E-01 +3.519379E-03 +2.306630E+00 +1.083425E+00 +7.695103E-02 +1.205473E-03 +1.872834E-01 +7.140475E-03 +7.077631E+00 +1.018246E+01 +8.321109E-02 +1.408010E-03 +2.025216E-01 +8.340388E-03 +2.094896E+01 +8.817186E+01 +3.233088E-02 +2.099749E-04 +8.000050E-02 +1.285635E-03 +1.356905E+01 +3.683229E+01 +1.859858E-01 +6.921278E-03 +5.173102E-01 +5.354619E-02 +5.056029E+00 +5.118575E+00 +1.910007E+00 +7.308066E-01 +4.648575E+00 +4.328846E+00 +2.856363E+00 +1.634636E+00 +4.328679E-01 +3.755784E-02 +1.053514E+00 +2.224694E-01 +1.692444E+00 +5.793158E-01 +5.559366E-02 +6.243434E-04 +1.353038E-01 +3.698224E-03 +2.368251E+00 +1.142832E+00 +7.949453E-02 +1.285913E-03 +1.934738E-01 +7.616955E-03 +7.118354E+00 +1.029753E+01 +8.426801E-02 +1.442436E-03 +2.050940E-01 +8.544309E-03 +2.046889E+01 +8.419814E+01 +3.182500E-02 +2.035334E-04 +7.874874E-02 +1.246195E-03 +1.326056E+01 +3.517079E+01 +1.833225E-01 +6.723237E-03 +5.099023E-01 +5.201406E-02 diff --git a/tests/regression_tests/random_ray_linear/test.py b/tests/regression_tests/random_ray_linear/test.py new file mode 100644 index 00000000000..10262789deb --- /dev/null +++ b/tests/regression_tests/random_ray_linear/test.py @@ -0,0 +1,26 @@ +import os + +import openmc +from openmc.examples import random_ray_lattice +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) + + +@pytest.mark.parametrize("shape", ["linear", "linear_xy"]) +def test_random_ray_source(shape): + with change_directory(shape): + openmc.reset_auto_ids() + model = random_ray_lattice() + model.settings.random_ray['source_shape'] = shape + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main()