diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 38a49b6021..d75a64d662 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -75,6 +75,7 @@ jobs: LIBMESH: ${{ matrix.libmesh }} NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX" OPENBLAS_NUM_THREADS: 1 + PYTEST_ADDOPTS: --cov=openmc --cov-report=lcov:coverage-python.lcov # libfabric complains about fork() as a result of using Python multiprocessing. # We can work around it with RDMAV_FORK_SAFE=1 in libfabric < 1.13 and with # FI_EFA_FORK_SAFE=1 in more recent versions. @@ -171,11 +172,37 @@ jobs: uses: mxschmitt/action-tmate@v3 timeout-minutes: 10 - - name: after_success + - name: Generate C++ coverage (gcovr) shell: bash run: | - cpp-coveralls -i src -i include -e src/external --exclude-pattern "/usr/*" --dump cpp_cov.json - coveralls --merge=cpp_cov.json --service=github + # Produce LCOV directly from gcov data in the build tree + gcovr \ + --root "$GITHUB_WORKSPACE" \ + --object-directory "$GITHUB_WORKSPACE/build" \ + --filter "$GITHUB_WORKSPACE/src" \ + --filter "$GITHUB_WORKSPACE/include" \ + --exclude "$GITHUB_WORKSPACE/src/external/.*" \ + --exclude "$GITHUB_WORKSPACE/src/include/openmc/external/.*" \ + --gcov-ignore-errors source_not_found \ + --gcov-ignore-errors output_error \ + --gcov-ignore-parse-errors suspicious_hits.warn \ + --print-summary \ + --lcov -o coverage-cpp.lcov || true + + - name: Merge C++ and Python coverage + shell: bash + run: | + # Merge C++ and Python LCOV into a single file for upload + cat coverage-cpp.lcov coverage-python.lcov > coverage.lcov + + - name: Upload coverage to Coveralls + if: ${{ hashFiles('coverage.lcov') != '' }} + uses: coverallsapp/github-action@v2 + with: + github-token: ${{ secrets.GITHUB_TOKEN }} + parallel: true + flag-name: C++ and Python + path-to-lcov: coverage.lcov finish: needs: main @@ -184,5 +211,5 @@ jobs: - name: Coveralls Finished uses: coverallsapp/github-action@v2 with: - github-token: ${{ secrets.github_token }} + github-token: ${{ secrets.GITHUB_TOKEN }} parallel-finished: true diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index 89c76525f8..74413e7b61 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -56,6 +56,27 @@ attributes: .. _io_chain_reaction: +-------------------- +```` Element +-------------------- + +The ```` element represents photon and electron sources associated with +the decay of a nuclide and contains information to construct an +:class:`openmc.stats.Univariate` object that represents this emission as an +energy distribution. This element has the following attributes: + + :type: + The type of :class:`openmc.stats.Univariate` source term. + + :particle: + The type of particle emitted, e.g., 'photon' or 'electron' + + :parameters: + The parameters of the source term, e.g., for a + :class:`openmc.stats.Discrete` source, the energies (in [eV]) at which the + particles are emitted and their relative intensities in [Bq/atom] (in other + words, decay constants). + ---------------------- ```` Element ---------------------- diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 26673faac2..720846c851 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -178,6 +178,16 @@ history-based parallelism. *Default*: false +-------------------------------- +```` Element +-------------------------------- + +The ```` element specifies the energy multiplier, expressed +in units of :math:`kT`, that determines when the free gas scattering approach is +used for elastic scattering. Values must be positive. + + *Default*: 400.0 + ----------------------------------- ```` Element ----------------------------------- diff --git a/docs/source/methods/charged_particles_physics.rst b/docs/source/methods/charged_particles_physics.rst new file mode 100644 index 0000000000..5d763074fd --- /dev/null +++ b/docs/source/methods/charged_particles_physics.rst @@ -0,0 +1,362 @@ +.. _methods_charged_particle_physics: + +======================== +Charged Particle Physics +======================== + +OpenMC neglects the spatial transport of charged particles (electrons and +positrons), assuming they deposit all their energy locally and produce +bremsstrahlung photons at their birth location. This approximation, called +thick-target bremsstrahlung (TTB) approximation is justified by the fact that +charged particles have much shorter stopping ranges compared to neutrons and +photons, especially in high-density materials. + +----------------------------- +Charged Particle Interactions +----------------------------- + +Bremsstrahlung +-------------- + +When a charged particle is decelerated in the field of an atom, some of its +kinetic energy is converted into electromagnetic radiation known as +bremsstrahlung, or 'braking radiation'. In each event, an electron or positron +with kinetic energy :math:`T` generates a photon with an energy :math:`E` +between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section +that is differential in photon energy, in the direction of the emitted photon, +and in the final direction of the charged particle. However, in Monte Carlo +simulations it is typical to integrate over the angular variables to obtain a +single differential cross section with respect to photon energy, which is often +expressed in the form + +.. math:: + :label: bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E} + \chi(Z, T, \kappa), + +where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, +\kappa)` is the scaled bremsstrahlung cross section, which is experimentally +measured. + +Because electrons are attracted to atomic nuclei whereas positrons are +repulsed, the cross section for positrons is smaller, though it approaches that +of electrons in the high energy limit. To obtain the positron cross section, we +multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used +in Salvat_, + +.. math:: + :label: positron-factor + + \begin{aligned} + F_{\text{p}}(Z,T) = + & 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\ + & + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\ + & - 1.8080\times 10^{-6}t^7), + \end{aligned} + +where + +.. math:: + :label: positron-factor-t + + t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right). + +:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for +positrons and electrons. Stopping power describes the average energy loss per +unit path length of a charged particle as it passes through matter: + +.. math:: + :label: stopping-power + + -\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T), + +where :math:`n` is the number density of the material and :math:`d\sigma/dE` is +the cross section differential in energy loss. The total stopping power +:math:`S(T)` can be separated into two components: the radiative stopping +power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to +bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`, +which refers to the energy loss due to inelastic collisions with bound +electrons in the material that result in ionization and excitation. The +radiative stopping power for electrons is given by + +.. math:: + :label: radiative-stopping-power + + S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa) + d\kappa. + + +To obtain the radiative stopping power for positrons, +:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`. + +While the models for photon interactions with matter described above can safely +assume interactions occur with free atoms, sampling the target atom based on +the macroscopic cross sections, molecular effects cannot necessarily be +disregarded for charged particle treatment. For compounds and mixtures, the +bremsstrahlung cross section is calculated using Bragg's additivity rule as + +.. math:: + :label: material-bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i + \chi(Z_i, T, \kappa), + +where the sum is over the constituent elements and :math:`\gamma_i` is the +atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping +power is calculated using Bragg's additivity rule as + +.. math:: + :label: material-radiative-stopping-power + + S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T), + +where :math:`w_i` is the mass fraction of the :math:`i`-th element and +:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using +:eq:`radiative-stopping-power`. The collision stopping power, however, is a +function of certain quantities such as the mean excitation energy :math:`I` and +the density effect correction :math:`\delta_F` that depend on molecular +properties. These quantities cannot simply be summed over constituent elements +in a compound, but should instead be calculated for the material. The Bethe +formula can be used to find the collision stopping power of the material: + +.. math:: + :label: material-collision-stopping-power + + S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M} + [\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)], + +where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass, +:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For +electrons, + +.. math:: + :label: F-electron + + F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2], + +while for positrons + +.. math:: + :label: F-positron + + F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 + + 4/(\tau + 2)^3]. + +The density effect correction :math:`\delta_F` takes into account the reduction +of the collision stopping power due to the polarization of the material the +charged particle is passing through by the electric field of the particle. +It can be evaluated using the method described by Sternheimer_, where the +equation for :math:`\delta_F` is + +.. math:: + :label: density-effect-correction + + \delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] - + l^2(1-\beta^2). + +Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition, +given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in +the :math:`i`-th subshell. The frequency :math:`l` is the solution of the +equation + +.. math:: + :label: density-effect-l + + \frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2}, + +where :math:`\bar{v}_i` is defined as + +.. math:: + :label: density-effect-nubar + + \bar{\nu}_i = h\nu_i \rho / h\nu_p. + +The plasma energy :math:`h\nu_p` of the medium is given by + +.. math:: + :label: plasma-frequency + + h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}}, + +where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the +material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator +energy, and :math:`\rho` is an adjustment factor introduced to give agreement +between the experimental values of the oscillator energies and the mean +excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are +defined as + +.. math:: + :label: density-effect-li + + \begin{aligned} + l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\ + l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0, + \end{aligned} + +where the second case applies to conduction electrons. For a conductor, +:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective +number of conduction electrons, and :math:`v_n = 0`. The adjustment factor +:math:`\rho` is determined using the equation for the mean excitation energy: + +.. math:: + :label: mean-excitation-energy + + \ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} + + f_n \ln (h\nu_pf_n^{1/2}). + +.. _ttb: + + +Thick-Target Bremsstrahlung Approximation ++++++++++++++++++++++++++++++++++++++++++ + +Since charged particles lose their energy on a much shorter distance scale than +neutral particles, not much error should be introduced by neglecting to +transport electrons. However, the bremsstrahlung emitted from high energy +electrons and positrons can travel far from the interaction site. Thus, even +without a full electron transport mode it is necessary to model bremsstrahlung. +We use a thick-target bremsstrahlung (TTB) approximation based on the models in +Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes +the charged particle loses all its energy in a single homogeneous material +region. + +To model bremsstrahlung using the TTB approximation, we need to know the number +of photons emitted by the charged particle and the energy distribution of the +photons. These quantities can be calculated using the continuous slowing down +approximation (CSDA). The CSDA assumes charged particles lose energy +continuously along their trajectory with a rate of energy loss equal to the +total stopping power, ignoring fluctuations in the energy loss. The +approximation is useful for expressing average quantities that describe how +charged particles slow down in matter. For example, the CSDA range approximates +the average path length a charged particle travels as it slows to rest: + +.. math:: + :label: csda-range + + R(T) = \int^T_0 \frac{dT'}{S(T')}. + +Actual path lengths will fluctuate around :math:`R(T)`. The average number of +photons emitted per unit path length is given by the inverse bremsstrahlung +mean free path: + +.. math:: + :label: inverse-bremsstrahlung-mfp + + \lambda_{\text{br}}^{-1}(T,E_{\text{cut}}) + = n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE + = n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa} + \chi(Z,T,\kappa)d\kappa. + +The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero +because the bremsstrahlung differential cross section diverges for small photon +energies but is finite for photon energies above some cutoff energy +:math:`E_{\text{cut}}`. The mean free path +:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the +photon number yield, defined as the average number of photons emitted with +energy greater than :math:`E_{\text{cut}}` as the charged particle slows down +from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is +given by + +.. math:: + :label: photon-number-yield + + Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})} + \lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T + \frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'. + +:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of +bremsstrahlung photons: the number of photons created with energy between +:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy +:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`. + +To simulate the emission of bremsstrahlung photons, the total stopping power +and bremsstrahlung differential cross section for positrons and electrons must +be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and +:eq:`material-radiative-stopping-power`. These quantities are used to build the +tabulated bremsstrahlung energy PDF and CDF for that material for each incident +energy :math:`T_k` on the energy grid. The following algorithm is then applied +to sample the photon energies: + +1. For an incident charged particle with energy :math:`T`, sample the number of + emitted photons as + + .. math:: + + N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor. + +2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1` + for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use + the composition method and sample from the PDF at either :math:`k` or + :math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can + be expressed as + + .. math:: + + p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1} + p_{\text{br}}(T_{k+1},E), + + where the interpolation weights are + + .. math:: + + \pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~ + \pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}. + + Sample either the index :math:`i = k` or :math:`i = k+1` according to the + point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`. + +3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`. + +3. Sample the photon energies using the inverse transform method with the + tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e., + + .. math:: + + E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} - + P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1 + \right]^{\frac{1}{1 + a_j}} + + where the interpolation factor :math:`a_j` is given by + + .. math:: + + a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)} + {\ln E_{j+1} - \ln E_j} + + and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le + P_{\text{br}}(T_i, E_{j+1})`. + +We ignore the range of the electron or positron, i.e., the bremsstrahlung +photons are produced in the same location that the charged particle was +created. The direction of the photons is assumed to be the same as the +direction of the incident charged particle, which is a reasonable approximation +at higher energies when the bremsstrahlung radiation is emitted at small +angles. + + +Electron-Positron Annihilation +------------------------------ + +When a positron collides with an electron, both particles are annihilated and +generally two photons with equal energy are created. If the kinetic energy of +the positron is high enough, the two photons can have different energies, and +the higher-energy photon is emitted preferentially in the direction of flight +of the positron. It is also possible to produce a single photon if the +interaction occurs with a bound electron, and in some cases three (or, rarely, +even more) photons can be emitted. However, the annihilation cross section is +largest for low-energy positrons, and as the positron energy decreases, the +angular distribution of the emitted photons becomes isotropic. + +In OpenMC, we assume the most likely case in which a low-energy positron (which +has already lost most of its energy to bremsstrahlung radiation) interacts with +an electron which is free and at rest. Two photons with energy equal to the +electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically +in opposite directions. + + +.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf + +.. _Salvat: https://doi.org/10.1787/32da5043-en + +.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067 diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 75c421c877..121d04b1de 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -14,6 +14,7 @@ Theory and Methodology random_numbers neutron_physics photon_physics + charged_particles_physics tallies eigenvalue depletion @@ -21,4 +22,4 @@ Theory and Methodology parallelization cmfd variance_reduction - random_ray \ No newline at end of file + random_ray diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 22d2c7f26a..d2bd3ac760 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -667,342 +667,6 @@ and Auger electrons: 5. Repeat from step 1 for vacancy left by the transition electron. -Electron-Positron Annihilation ------------------------------- - -When a positron collides with an electron, both particles are annihilated and -generally two photons with equal energy are created. If the kinetic energy of -the positron is high enough, the two photons can have different energies, and -the higher-energy photon is emitted preferentially in the direction of flight -of the positron. It is also possible to produce a single photon if the -interaction occurs with a bound electron, and in some cases three (or, rarely, -even more) photons can be emitted. However, the annihilation cross section is -largest for low-energy positrons, and as the positron energy decreases, the -angular distribution of the emitted photons becomes isotropic. - -In OpenMC, we assume the most likely case in which a low-energy positron (which -has already lost most of its energy to bremsstrahlung radiation) interacts with -an electron which is free and at rest. Two photons with energy equal to the -electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically -in opposite directions. - -Bremsstrahlung --------------- - -When a charged particle is decelerated in the field of an atom, some of its -kinetic energy is converted into electromagnetic radiation known as -bremsstrahlung, or 'braking radiation'. In each event, an electron or positron -with kinetic energy :math:`T` generates a photon with an energy :math:`E` -between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section -that is differential in photon energy, in the direction of the emitted photon, -and in the final direction of the charged particle. However, in Monte Carlo -simulations it is typical to integrate over the angular variables to obtain a -single differential cross section with respect to photon energy, which is often -expressed in the form - -.. math:: - :label: bremsstrahlung-dcs - - \frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E} - \chi(Z, T, \kappa), - -where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, -\kappa)` is the scaled bremsstrahlung cross section, which is experimentally -measured. - -Because electrons are attracted to atomic nuclei whereas positrons are -repulsed, the cross section for positrons is smaller, though it approaches that -of electrons in the high energy limit. To obtain the positron cross section, we -multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used -in Salvat_, - -.. math:: - :label: positron-factor - - \begin{aligned} - F_{\text{p}}(Z,T) = - & 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\ - & + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\ - & - 1.8080\times 10^{-6}t^7), - \end{aligned} - -where - -.. math:: - :label: positron-factor-t - - t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right). - -:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for -positrons and electrons. Stopping power describes the average energy loss per -unit path length of a charged particle as it passes through matter: - -.. math:: - :label: stopping-power - - -\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T), - -where :math:`n` is the number density of the material and :math:`d\sigma/dE` is -the cross section differential in energy loss. The total stopping power -:math:`S(T)` can be separated into two components: the radiative stopping -power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to -bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`, -which refers to the energy loss due to inelastic collisions with bound -electrons in the material that result in ionization and excitation. The -radiative stopping power for electrons is given by - -.. math:: - :label: radiative-stopping-power - - S_{\text{rad}}(T) = n \frac{Z^2}{\beta^2} T \int_0^1 \chi(Z,T,\kappa) - d\kappa. - - -To obtain the radiative stopping power for positrons, -:eq:`radiative-stopping-power` is multiplied by :eq:`positron-factor`. - -While the models for photon interactions with matter described above can safely -assume interactions occur with free atoms, sampling the target atom based on -the macroscopic cross sections, molecular effects cannot necessarily be -disregarded for charged particle treatment. For compounds and mixtures, the -bremsstrahlung cross section is calculated using Bragg's additivity rule as - -.. math:: - :label: material-bremsstrahlung-dcs - - \frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i - \chi(Z_i, T, \kappa), - -where the sum is over the constituent elements and :math:`\gamma_i` is the -atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping -power is calculated using Bragg's additivity rule as - -.. math:: - :label: material-radiative-stopping-power - - S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T), - -where :math:`w_i` is the mass fraction of the :math:`i`-th element and -:math:`S_{\text{rad},i}(T)` is found for element :math:`i` using -:eq:`radiative-stopping-power`. The collision stopping power, however, is a -function of certain quantities such as the mean excitation energy :math:`I` and -the density effect correction :math:`\delta_F` that depend on molecular -properties. These quantities cannot simply be summed over constituent elements -in a compound, but should instead be calculated for the material. The Bethe -formula can be used to find the collision stopping power of the material: - -.. math:: - :label: material-collision-stopping-power - - S_{\text{col}}(T) = \frac{2 \pi r_e^2 m_e c^2}{\beta^2} N_A \frac{Z}{A_M} - [\ln(T^2/I^2) + \ln(1 + \tau/2) + F(\tau) - \delta_F(T)], - -where :math:`N_A` is Avogadro's number, :math:`A_M` is the molar mass, -:math:`\tau = T/m_e`, and :math:`F(\tau)` depends on the particle type. For -electrons, - -.. math:: - :label: F-electron - - F_{-}(\tau) = (1 - \beta^2)[1 + \tau^2/8 - (2\tau + 1) \ln2], - -while for positrons - -.. math:: - :label: F-positron - - F_{+}(\tau) = 2\ln2 - (\beta^2/12)[23 + 14/(\tau + 2) + 10/(\tau + 2)^2 + - 4/(\tau + 2)^3]. - -The density effect correction :math:`\delta_F` takes into account the reduction -of the collision stopping power due to the polarization of the material the -charged particle is passing through by the electric field of the particle. -It can be evaluated using the method described by Sternheimer_, where the -equation for :math:`\delta_F` is - -.. math:: - :label: density-effect-correction - - \delta_F(\beta) = \sum_{i=1}^n f_i \ln[(l_i^2 + l^2)/l_i^2] - - l^2(1-\beta^2). - -Here, :math:`f_i` is the oscillator strength of the :math:`i`-th transition, -given by :math:`f_i = n_i/Z`, where :math:`n_i` is the number of electrons in -the :math:`i`-th subshell. The frequency :math:`l` is the solution of the -equation - -.. math:: - :label: density-effect-l - - \frac{1}{\beta^2} - 1 = \sum_{i=1}^{n} \frac{f_i}{\bar{\nu}_i^2 + l^2}, - -where :math:`\bar{v}_i` is defined as - -.. math:: - :label: density-effect-nubar - - \bar{\nu}_i = h\nu_i \rho / h\nu_p. - -The plasma energy :math:`h\nu_p` of the medium is given by - -.. math:: - :label: plasma-frequency - - h\nu_p = \sqrt{\frac{(hc)^2 r_e \rho_m N_A Z}{\pi A}}, - -where :math:`A` is the atomic weight and :math:`\rho_m` is the density of the -material. In :eq:`density-effect-nubar`, :math:`h\nu_i` is the oscillator -energy, and :math:`\rho` is an adjustment factor introduced to give agreement -between the experimental values of the oscillator energies and the mean -excitation energy. The :math:`l_i` in :eq:`density-effect-correction` are -defined as - -.. math:: - :label: density-effect-li - - \begin{aligned} - l_i &= (\bar{\nu}_i^2 + 2/3f_i)^{1/2} ~~~~&\text{for}~~ \bar{\nu}_i > 0 \\ - l_n &= f_n^{1/2} ~~~~&\text{for}~~ \bar{\nu}_n = 0, - \end{aligned} - -where the second case applies to conduction electrons. For a conductor, -:math:`f_n` is given by :math:`n_c/Z`, where :math:`n_c` is the effective -number of conduction electrons, and :math:`v_n = 0`. The adjustment factor -:math:`\rho` is determined using the equation for the mean excitation energy: - -.. math:: - :label: mean-excitation-energy - - \ln I = \sum_{i=1}^{n-1} f_i \ln[(h\nu_i\rho)^2 + 2/3f_i(h\nu_p)^2]^{1/2} + - f_n \ln (h\nu_pf_n^{1/2}). - -.. _ttb: - -Thick-Target Bremsstrahlung Approximation -+++++++++++++++++++++++++++++++++++++++++ - -Since charged particles lose their energy on a much shorter distance scale than -neutral particles, not much error should be introduced by neglecting to -transport electrons. However, the bremsstrahlung emitted from high energy -electrons and positrons can travel far from the interaction site. Thus, even -without a full electron transport mode it is necessary to model bremsstrahlung. -We use a thick-target bremsstrahlung (TTB) approximation based on the models in -Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes -the charged particle loses all its energy in a single homogeneous material -region. - -To model bremsstrahlung using the TTB approximation, we need to know the number -of photons emitted by the charged particle and the energy distribution of the -photons. These quantities can be calculated using the continuous slowing down -approximation (CSDA). The CSDA assumes charged particles lose energy -continuously along their trajectory with a rate of energy loss equal to the -total stopping power, ignoring fluctuations in the energy loss. The -approximation is useful for expressing average quantities that describe how -charged particles slow down in matter. For example, the CSDA range approximates -the average path length a charged particle travels as it slows to rest: - -.. math:: - :label: csda-range - - R(T) = \int^T_0 \frac{dT'}{S(T')}. - -Actual path lengths will fluctuate around :math:`R(T)`. The average number of -photons emitted per unit path length is given by the inverse bremsstrahlung -mean free path: - -.. math:: - :label: inverse-bremsstrahlung-mfp - - \lambda_{\text{br}}^{-1}(T,E_{\text{cut}}) - = n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE - = n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa} - \chi(Z,T,\kappa)d\kappa. - -The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero -because the bremsstrahlung differential cross section diverges for small photon -energies but is finite for photon energies above some cutoff energy -:math:`E_{\text{cut}}`. The mean free path -:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the -photon number yield, defined as the average number of photons emitted with -energy greater than :math:`E_{\text{cut}}` as the charged particle slows down -from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is -given by - -.. math:: - :label: photon-number-yield - - Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})} - \lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T - \frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'. - -:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of -bremsstrahlung photons: the number of photons created with energy between -:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy -:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`. - -To simulate the emission of bremsstrahlung photons, the total stopping power -and bremsstrahlung differential cross section for positrons and electrons must -be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and -:eq:`material-radiative-stopping-power`. These quantities are used to build the -tabulated bremsstrahlung energy PDF and CDF for that material for each incident -energy :math:`T_k` on the energy grid. The following algorithm is then applied -to sample the photon energies: - -1. For an incident charged particle with energy :math:`T`, sample the number of - emitted photons as - - .. math:: - - N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor. - -2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1` - for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use - the composition method and sample from the PDF at either :math:`k` or - :math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can - be expressed as - - .. math:: - - p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1} - p_{\text{br}}(T_{k+1},E), - - where the interpolation weights are - - .. math:: - - \pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~ - \pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}. - - Sample either the index :math:`i = k` or :math:`i = k+1` according to the - point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`. - -3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`. - -3. Sample the photon energies using the inverse transform method with the - tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e., - - .. math:: - - E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} - - P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1 - \right]^{\frac{1}{1 + a_j}} - - where the interpolation factor :math:`a_j` is given by - - .. math:: - - a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)} - {\ln E_{j+1} - \ln E_j} - - and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le - P_{\text{br}}(T_i, E_{j+1})`. - -We ignore the range of the electron or positron, i.e., the bremsstrahlung -photons are produced in the same location that the charged particle was -created. The direction of the photons is assumed to be the same as the -direction of the incident charged particle, which is a reasonable approximation -at higher energies when the bremsstrahlung radiation is emitted at small -angles. .. _photon_production: @@ -1070,5 +734,3 @@ emitted photon. .. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf .. _Salvat: https://doi.org/10.1787/32da5043-en - -.. _Sternheimer: https://doi.org/10.1103/PhysRevB.26.6067 diff --git a/include/openmc/constants.h b/include/openmc/constants.h index df13da3707..a0d1646131 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5}; // prevent extremely large adjoint source terms from being generated. constexpr double ZERO_FLUX_CUTOFF {1e-22}; +// The minimum macroscopic cross section value considered non-void for the +// random ray solver. Materials with any group with a cross section below this +// value will be converted to pure void. +constexpr double MINIMUM_MACRO_XS {1e-6}; + // ============================================================================ // MATH AND PHYSICAL CONSTANTS diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 329c776d03..60b88a153b 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -164,8 +164,8 @@ namespace data { // Minimum/maximum transport energy for each particle type. Order corresponds to // that of the ParticleType enum -extern array energy_min; -extern array energy_max; +extern array energy_min; +extern array energy_max; //! Minimum temperature in [K] that nuclide data is available at extern double temperature_min; diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 1df2b3ef19..e8188e093a 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -154,9 +154,10 @@ struct NuclideMicroXS { // Energy and temperature last used to evaluate these cross sections. If // these values have changed, then the cross sections must be re-evaluated. - double last_E {0.0}; //!< Last evaluated energy - double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant - //!< * temperature (eV)) + double last_E {0.0}; //!< Last evaluated energy + double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant + //!< * temperature (eV)) + double ncrystal_xs {-1.0}; //!< NCrystal cross section }; //============================================================================== diff --git a/include/openmc/physics.h b/include/openmc/physics.h index f62f43a02f..2472d97993 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -10,13 +10,6 @@ namespace openmc { -//============================================================================== -// Constants -//============================================================================== - -// Monoatomic ideal-gas scattering treatment threshold -constexpr double FREE_GAS_THRESHOLD {400.0}; - //============================================================================== // Non-member functions //============================================================================== diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 78351fcc5f..d4e8027346 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -27,8 +27,9 @@ public: //---------------------------------------------------------------------------- // Methods - virtual void update_neutron_source(double k_eff); - double compute_k_eff(double k_eff_old) const; + virtual void update_single_neutron_source(SourceRegionHandle& srh); + virtual void update_all_neutron_sources(); + void compute_k_eff(); virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); @@ -41,7 +42,7 @@ public: void output_to_vtk() const; void convert_external_sources(); void count_external_source_regions(); - void set_adjoint_sources(const vector& forward_flux); + void set_adjoint_sources(); void flux_swap(); virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; @@ -54,9 +55,8 @@ public: bool is_target_void); void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx, int32_t target_material_id, bool is_target_void); - void prepare_base_source_regions(); SourceRegionHandle get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u); + SourceRegionKey sr_key, Position r, Direction u); void finalize_discovered_source_regions(); void apply_transport_stabilization(); int64_t n_source_regions() const @@ -67,6 +67,10 @@ public: { return source_regions_.n_source_regions() * negroups_; } + int64_t lookup_base_source_region_idx(const GeometryState& p) const; + SourceRegionKey lookup_source_region_key(const GeometryState& p) const; + int64_t lookup_mesh_bin(int64_t sr, Position r) const; + int lookup_mesh_idx(int64_t sr) const; //---------------------------------------------------------------------------- // Static Data members @@ -86,6 +90,7 @@ public: //---------------------------------------------------------------------------- // Public Data members + double k_eff_ {1.0}; // Eigenvalue bool mapped_all_tallies_ {false}; // If all source regions have been visited int64_t n_external_source_regions_ {0}; // Total number of source regions with @@ -110,14 +115,6 @@ public: // The abstract container holding all source region-specific data SourceRegionContainer source_regions_; - // Base source region container. When source region subdivision via mesh - // is in use, this container holds the original (non-subdivided) material - // filled cell instance source regions. These are useful as they can be - // initialized with external source and mesh domain information ahead of time. - // Then, dynamically discovered source regions can be initialized by cloning - // their base region. - SourceRegionContainer base_source_regions_; - // Parallel hash map holding all source regions discovered during // a single iteration. This is a threadsafe data structure that is cleaned // out after each iteration and stored in the "source_regions_" container. @@ -134,8 +131,17 @@ public: // Map that relates a SourceRegionKey to the external source index. This map // is used to check if there are any point sources within a subdivided source // region at the time it is discovered. - std::unordered_map - point_source_map_; + std::unordered_map, SourceRegionKey::HashFunctor> + external_point_source_map_; + + // Map that relates a base source region index to the external source index. + // This map is used to check if there are any volumetric sources within a + // subdivided source region at the time it is discovered. + std::unordered_map> external_volumetric_source_map_; + + // Map that relates a base source region index to a mesh index. This map + // is used to check which subdivision mesh is present in a source region. + std::unordered_map mesh_map_; // If transport corrected MGXS data is being used, there may be negative // in-group scattering cross sections that can result in instability in MOC @@ -147,12 +153,11 @@ protected: //---------------------------------------------------------------------------- // Methods void apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh); - void apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, - const vector& instances); - void apply_external_source_to_cell_and_children(int32_t i_cell, - Discrete* discrete, double strength_factor, int32_t target_material_id); + int src_idx, SourceRegionHandle& srh); + void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx, + int target_material_id, const vector& instances); + void apply_external_source_to_cell_and_children( + int32_t i_cell, int src_idx, int32_t target_material_id); virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g); void set_flux_to_source(int64_t sr, int g); virtual void set_flux_to_old_flux(int64_t sr, int g); diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 67fdd99f88..0098c78200 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain { public: //---------------------------------------------------------------------------- // Methods - void update_neutron_source(double k_eff) override; + void update_single_neutron_source(SourceRegionHandle& srh) override; void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration) override; diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index abf2a26881..40c67ef954 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -48,7 +48,6 @@ public: 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 - static bool mesh_subdivision_enabled_; // Flag for mesh subdivision static RandomRaySampleMethod sample_method_; // Flag for sampling method //---------------------------------------------------------------------------- diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b94e7401b3..3dec48bf26 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -21,11 +21,7 @@ public: // Methods void compute_segment_correction_factors(); void apply_fixed_sources_and_mesh_domains(); - void prepare_fixed_sources_adjoint(vector& forward_flux, - SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, - std::unordered_map& - forward_source_region_map); + void prepare_fixed_sources_adjoint(); void simulate(); void output_simulation_results() const; void instability_check( @@ -45,9 +41,6 @@ private: // Contains all flat source region data unique_ptr domain_; - // Random ray eigenvalue - double k_eff_ {1.0}; - // Tracks the average FSR miss rate for analysis and reporting double avg_miss_rate_ {0.0}; diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 5c5b31f392..0f5a747fff 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -308,7 +308,6 @@ public: //---------------------------------------------------------------------------- // Constructors SourceRegion(int negroups, bool is_linear); - SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr); SourceRegion() = default; //---------------------------------------------------------------------------- diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 999b1b83f8..78bfa088e6 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -147,6 +147,8 @@ extern std::unordered_set source_write_surf_id; //!< Surface ids where sources will be written extern double source_rejection_fraction; //!< Minimum fraction of source sites //!< that must be accepted +extern double free_gas_threshold; //!< Threshold multiplier for free gas + //!< scattering treatment extern int max_history_splits; //!< maximum number of particle splits for weight windows diff --git a/openmc/data/data.py b/openmc/data/data.py index 2142a5dc90..5ecadd37be 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -324,7 +324,7 @@ def atomic_mass(isotope): # isotopes of their element (e.g. C0), calculate the atomic mass as # the sum of the atomic mass times the natural abundance of the isotopes # that make up the element. - for element in ['C', 'Zn', 'Pt', 'Os', 'Tl']: + for element in ['C', 'Zn', 'Pt', 'Os', 'Tl', 'V']: isotope_zero = element.lower() + '0' _ATOMIC_MASS[isotope_zero] = 0. for iso, abundance in isotopes(element): diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 1a11d3614f..c8a0bb5e7e 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -591,7 +591,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: openmc.stats.Univariate or None Distribution of energies in [eV] of photons emitted from decay, or None if no photon source exists. Note that the probabilities represent - intensities, given as [Bq]. + intensities, given as [Bq/atom] (in other words, decay constants). """ if not _DECAY_PHOTON_ENERGY: chain_file = openmc.config.get('chain_file') diff --git a/openmc/deplete/d1s.py b/openmc/deplete/d1s.py index 311bc2d69c..f51dea4160 100644 --- a/openmc/deplete/d1s.py +++ b/openmc/deplete/d1s.py @@ -108,14 +108,15 @@ def time_correction_factors( # Create a 2D array for the time correction factors h = np.zeros((n_timesteps, n_nuclides)) - for i, (dt, rate) in enumerate(zip(timesteps, source_rates)): - # Precompute the exponential terms. Since (1 - exp(-x)) is susceptible to - # roundoff error, use expm1 instead (which computes exp(x) - 1) - g = np.exp(-decay_rate*dt) - one_minus_g = -np.expm1(-decay_rate*dt) + # Precompute all exponential terms with same shape as h + decay_dt = decay_rate[np.newaxis, :] * timesteps[:, np.newaxis] + g = np.exp(-decay_dt) + one_minus_g = -np.expm1(-decay_dt) + # Apply recurrence relation step by step + for i in range(len(timesteps)): # Eq. (4) in doi:10.1016/j.fusengdes.2019.111399 - h[i + 1] = rate*one_minus_g + h[i]*g + h[i + 1] = source_rates[i] * one_minus_g[i] + h[i] * g[i] return {nuclides[i]: h[:, i] for i in range(n_nuclides)} diff --git a/openmc/plotter.py b/openmc/plotter.py index 693cdaca4e..abd8ab6dd4 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -501,7 +501,7 @@ def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None, elif ncrystal_cfg: import NCrystal nc_scatter = NCrystal.createScatter(ncrystal_cfg) - nc_func = nc_scatter.crossSectionNonOriented + nc_func = nc_scatter.xsect nc_emax = 5 # eV # this should be obtained from NCRYSTAL_MAX_ENERGY energy_grid = np.union1d(np.geomspace(min(energy_grid), 1.1*nc_emax, diff --git a/openmc/settings.py b/openmc/settings.py index 7bbc31c192..2f8a2b1248 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -84,6 +84,10 @@ class Settings: history-based parallelism. .. versionadded:: 0.12 + free_gas_threshold : float + Energy multiplier (in units of :math:`kT`) below which the free gas + scattering treatment is applied for elastic scattering. If not + specified, a value of 400.0 is used. generations_per_batch : int Number of generations per batch ifp_n_generation : int @@ -376,6 +380,7 @@ class Settings: self._seed = None self._stride = None self._survival_biasing = None + self._free_gas_threshold = None # Shannon entropy mesh self._entropy_mesh = None @@ -1255,6 +1260,17 @@ class Settings: cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1) self._source_rejection_fraction = source_rejection_fraction + @property + def free_gas_threshold(self) -> float | None: + return self._free_gas_threshold + + @free_gas_threshold.setter + def free_gas_threshold(self, free_gas_threshold: float | None): + if free_gas_threshold is not None: + cv.check_type('free gas threshold', free_gas_threshold, Real) + cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0) + self._free_gas_threshold = free_gas_threshold + def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1641,6 +1657,7 @@ class Settings: if mesh_memo is not None: mesh_memo.add(ww.mesh.id) + def _create_weight_windows_on_subelement(self, root): if self._weight_windows_on is not None: elem = ET.SubElement(root, "weight_windows_on") elem.text = str(self._weight_windows_on).lower() @@ -1714,9 +1731,15 @@ class Settings: 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: + # See if a element already exists -- if not, add it + path = f"./mesh[@id='{mesh.id}']" + if root.find(path) is None: root.append(mesh.to_xml_element()) - mesh_memo.add(mesh.id) + if mesh_memo is not None: + mesh_memo.add(mesh.id) + elif isinstance(value, bool): + subelement = ET.SubElement(element, key) + subelement.text = str(value).lower() else: subelement = ET.SubElement(element, key) subelement.text = str(value) @@ -1726,6 +1749,11 @@ class Settings: element = ET.SubElement(root, "source_rejection_fraction") element.text = str(self._source_rejection_fraction) + def _create_free_gas_threshold_subelement(self, root): + if self._free_gas_threshold is not None: + element = ET.SubElement(root, "free_gas_threshold") + element.text = str(self._free_gas_threshold) + def _eigenvalue_from_xml_element(self, root): elem = root.find('eigenvalue') if elem is not None: @@ -2074,10 +2102,16 @@ class Settings: ww = WeightWindows.from_xml_element(elem, meshes) self.weight_windows.append(ww) + def _weight_windows_on_from_xml_element(self, root): text = get_text(root, 'weight_windows_on') if text is not None: self.weight_windows_on = text in ('true', '1') + def _weight_windows_file_from_xml_element(self, root): + text = get_text(root, 'weight_windows_file') + if text is not None: + self.weight_windows_file = text + def _weight_window_checkpoints_from_xml_element(self, root): elem = root.find('weight_window_checkpoints') if elem is None: @@ -2103,7 +2137,7 @@ class Settings: if text is not None: self.max_tracks = int(text) - def _random_ray_from_xml_element(self, root): + def _random_ray_from_xml_element(self, root, meshes=None): elem = root.find('random_ray') if elem is not None: self.random_ray = {} @@ -2130,7 +2164,11 @@ class Settings: 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) + mesh_id = int(get_text(mesh_elem, 'id')) + if meshes and mesh_id in meshes: + mesh = meshes[mesh_id] + else: + mesh = MeshBase.from_xml_element(mesh_elem) domains = [] for domain_elem in mesh_elem.findall('domain'): domain_id = int(get_text(domain_elem, "id")) @@ -2154,6 +2192,11 @@ class Settings: if text is not None: self.source_rejection_fraction = float(text) + def _free_gas_threshold_from_xml_element(self, root): + text = get_text(root, 'free_gas_threshold') + if text is not None: + self.free_gas_threshold = float(text) + def to_xml_element(self, mesh_memo=None): """Create a 'settings' element to be written to an XML file. @@ -2214,6 +2257,7 @@ class Settings: self._create_log_grid_bins_subelement(element) self._create_write_initial_source_subelement(element) self._create_weight_windows_subelement(element, mesh_memo) + self._create_weight_windows_on_subelement(element) self._create_weight_window_generators_subelement(element, mesh_memo) self._create_weight_windows_file_element(element) self._create_weight_window_checkpoints_subelement(element) @@ -2223,6 +2267,7 @@ class Settings: self._create_random_ray_subelement(element, mesh_memo) self._create_use_decay_photons_subelement(element) self._create_source_rejection_fraction_subelement(element) + self._create_free_gas_threshold_subelement(element) # Clean the indentation in the file to be user-readable clean_indentation(element) @@ -2324,14 +2369,17 @@ class Settings: settings._log_grid_bins_from_xml_element(elem) settings._write_initial_source_from_xml_element(elem) settings._weight_windows_from_xml_element(elem, meshes) + settings._weight_windows_on_from_xml_element(elem) + settings._weight_windows_file_from_xml_element(elem) settings._weight_window_generators_from_xml_element(elem, meshes) settings._weight_window_checkpoints_from_xml_element(elem) settings._max_history_splits_from_xml_element(elem) settings._max_tracks_from_xml_element(elem) settings._max_secondaries_from_xml_element(elem) - settings._random_ray_from_xml_element(elem) + settings._random_ray_from_xml_element(elem, meshes) settings._use_decay_photons_from_xml_element(elem) settings._source_rejection_fraction_from_xml_element(elem) + settings._free_gas_threshold_from_xml_element(elem) return settings diff --git a/openmc/source.py b/openmc/source.py index c463ccb275..9b730cf1de 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -260,7 +260,7 @@ class IndependentSource(SourceBase): time distribution of source sites strength : float Strength of the source - particle : {'neutron', 'photon'} + particle : {'neutron', 'photon', 'electron', 'positron'} Source particle type domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe Domains to reject based on, i.e., if a sampled spatial location is not @@ -299,7 +299,7 @@ class IndependentSource(SourceBase): .. versionadded:: 0.14.0 - particle : {'neutron', 'photon'} + particle : {'neutron', 'photon', 'electron', 'positron'} Source particle type constraints : dict Constraints on sampled source particles. Valid keys include @@ -404,7 +404,8 @@ class IndependentSource(SourceBase): @particle.setter def particle(self, particle): - cv.check_value('source particle', particle, ['neutron', 'photon']) + cv.check_value('source particle', particle, + ['neutron', 'photon', 'electron', 'positron']) self._particle = particle def populate_xml_element(self, element): diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 29c11921cb..a763db3971 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -536,7 +536,7 @@ class StatePoint: def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None, exact_filters=False, exact_nuclides=False, exact_scores=False, - multiply_density=None, derivative=None): + multiply_density=None, derivative=None, filter_type=None): """Finds and returns a Tally object with certain properties. This routine searches the list of Tallies and returns the first Tally @@ -580,6 +580,9 @@ class StatePoint: to the same value as this parameter. derivative : openmc.TallyDerivative, optional TallyDerivative object to match. + filter_type : type, optional + If not None, the Tally must have at least one Filter that is an + instance of this type. For example `openmc.MeshFilter`. Returns ------- @@ -653,6 +656,10 @@ class StatePoint: if not contains_filters: continue + if filter_type is not None: + if not any(isinstance(f, filter_type) for f in test_tally.filters): + continue + # Determine if Tally has the queried Nuclide(s) if nuclides: if not all(nuclide in test_tally.nuclides for nuclide in nuclides): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index d6cf19f2b8..c48cc00757 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -295,6 +295,20 @@ class Discrete(Univariate): """ return np.sum(self.p) + def mean(self) -> float: + """Return mean of the discrete distribution + + The mean is the weighted average of the discrete values. + + .. versionadded:: 0.15.3 + + Returns + ------- + float + Mean of discrete distribution + """ + return np.sum(self.x * self.p) / np.sum(self.p) + def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Discrete: r"""Remove low-importance points from discrete distribution. @@ -413,6 +427,18 @@ class Uniform(Univariate): rng = np.random.RandomState(seed) return rng.uniform(self.a, self.b, n_samples) + def mean(self) -> float: + """Return mean of the uniform distribution + + .. versionadded:: 0.15.3 + + Returns + ------- + float + Mean of uniform distribution + """ + return 0.5 * (self.a + self.b) + def to_xml_element(self, element_name: str): """Return XML representation of the uniform distribution @@ -1123,7 +1149,7 @@ class Tabular(Univariate): """ interpolation = get_text(elem, 'interpolation') - params = get_elem_list(elem, "parameters", float) + params = get_elem_list(elem, "parameters", float) m = (len(params) + 1)//2 # +1 for when len(params) is odd x = params[:m] p = params[m:] @@ -1347,6 +1373,30 @@ class Mixture(Univariate): for p, dist in zip(self.probability, self.distribution) ]) + def mean(self) -> float: + """Return mean of the mixture distribution + + The mean is the weighted average of the means of the component + distributions, weighted by probability * integral. + + .. versionadded:: 0.15.3 + + Returns + ------- + float + Mean of the mixture distribution + """ + # Weight each component by its probability and integral + weights = [p*dist.integral() for p, dist in + zip(self.probability, self.distribution)] + total_weight = sum(weights) + + if total_weight == 0: + return 0.0 + + return sum([w*dist.mean() for w, dist in + zip(weights, self.distribution)]) / total_weight + def clip(self, tolerance: float = 1e-6, inplace: bool = False) -> Mixture: r"""Remove low-importance points / distributions @@ -1369,14 +1419,14 @@ class Mixture(Univariate): Distribution with low-importance points / distributions removed """ - # Determine integral of original distribution to compare later - original_integral = self.integral() + # Calculate mean * integral for original distribution to compare later. + original_mean_integral = self.mean() * self.integral() # Determine indices for any distributions that contribute non-negligibly - # to overall intensity - intensities = [prob*dist.integral() for prob, dist in - zip(self.probability, self.distribution)] - indices = _intensity_clip(intensities, tolerance=tolerance) + # to overall mean * integral + mean_integrals = [prob*dist.mean()*dist.integral() for prob, dist in + zip(self.probability, self.distribution)] + indices = _intensity_clip(mean_integrals, tolerance=tolerance) # Clip mixture of distributions probability = self.probability[indices] @@ -1397,12 +1447,14 @@ class Mixture(Univariate): # Create new distribution new_dist = type(self)(probability, distribution) - # Show warning if integral of new distribution is not within - # tolerance of original - diff = (original_integral - new_dist.integral())/original_integral + # Show warning if mean * integral of new distribution is not within + # tolerance of original. For energy distributions, mean * integral + # represents total energy. + new_mean_integral = new_dist.mean() * new_dist.integral() + diff = (original_mean_integral - new_mean_integral)/original_mean_integral if diff > tolerance: - warn("Clipping mixture distribution resulted in an integral that is " - f"lower by a fraction of {diff} when tolerance={tolerance}.") + warn("Clipping mixture distribution resulted in a mean*integral " + f"that is lower by a fraction of {diff} when tolerance={tolerance}.") return new_dist diff --git a/pyproject.toml b/pyproject.toml index 7705de8dae..2d67e83401 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,8 +48,15 @@ docs = [ "sphinxcontrib-svg2pdfconverter", "sphinx-rtd-theme" ] -test = ["packaging", "pytest", "pytest-cov", "colorama", "openpyxl"] -ci = ["cpp-coveralls", "coveralls"] +test = [ + "packaging", + "pytest", + "pytest-cov>=4.0", + "pytest-rerunfailures", + "colorama", + "openpyxl", +] +ci = ["coverage>=7.4", "gcovr>=7.2"] vtk = ["vtk"] [project.urls] diff --git a/src/finalize.cpp b/src/finalize.cpp index 25b471a8bd..659f390b39 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -85,6 +85,7 @@ int openmc_finalize() settings::time_cutoff = {INFTY, INFTY, INFTY, INFTY}; settings::entropy_on = false; settings::event_based = false; + settings::free_gas_threshold = 400.0; settings::gen_per_batch = 1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; @@ -155,8 +156,8 @@ int openmc_finalize() simulation::entropy_mesh = nullptr; simulation::ufs_mesh = nullptr; - data::energy_max = {INFTY, INFTY}; - data::energy_min = {0.0, 0.0}; + data::energy_max = {INFTY, INFTY, INFTY, INFTY}; + data::energy_min = {0.0, 0.0, 0.0, 0.0}; data::temperature_min = 0.0; data::temperature_max = INFTY; model::root_universe = -1; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 7cb84640d6..5ae6e30ee2 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -31,8 +31,8 @@ namespace openmc { //============================================================================== namespace data { -array energy_min {0.0, 0.0}; -array energy_max {INFTY, INFTY}; +array energy_min {0.0, 0.0, 0.0, 0.0}; +array energy_max {INFTY, INFTY, INFTY, INFTY}; double temperature_min {INFTY}; double temperature_max {0.0}; std::unordered_map nuclide_map; diff --git a/src/particle.cpp b/src/particle.cpp index 31b469adc1..748c698fc3 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -229,7 +229,7 @@ void Particle::event_advance() { // Sample a distance to collision if (type() == ParticleType::electron || type() == ParticleType::positron) { - collision_distance() = 0.0; + collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0; } else if (macro_xs().total == 0.0) { collision_distance() = INFINITY; } else { @@ -861,10 +861,12 @@ void Particle::update_neutron_xs( // If the cache doesn't match, recalculate micro xs if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT || - i_sab != micro.index_sab || sab_frac != micro.sab_frac) { + i_sab != micro.index_sab || sab_frac != micro.sab_frac || + ncrystal_xs != micro.ncrystal_xs) { data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this); // If NCrystal is being used, update micro cross section cache + micro.ncrystal_xs = ncrystal_xs; if (ncrystal_xs >= 0.0) { data::nuclides[i_nuclide]->calculate_elastic_xs(*this); ncrystal_update_micro(ncrystal_xs, micro); diff --git a/src/physics.cpp b/src/physics.cpp index e947fecbb9..3a17077b36 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -851,7 +851,7 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, // otherwise, use free gas model } else { - if (E >= FREE_GAS_THRESHOLD * kT && nuc.awr_ > 1.0) { + if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) { return {}; } else { sampling_method = ResScatMethod::cxs; diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 4092388308..1bf27e1eda 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -53,24 +53,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; source_regions_ = SourceRegionContainer(negroups_, is_linear); - source_regions_.assign( - base_source_regions, SourceRegion(negroups_, is_linear)); - - // Initialize materials - int64_t source_region_id = 0; - for (int i = 0; i < model::cells.size(); i++) { - Cell& cell = *model::cells[i]; - if (cell.type_ == Fill::MATERIAL) { - for (int j = 0; j < cell.n_instances(); j++) { - source_regions_.material(source_region_id++) = cell.material(j); - } - } - } - - // Sanity check - if (source_region_id != base_source_regions) { - fatal_error("Unexpected number of source regions"); - } // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -118,34 +100,24 @@ void FlatSourceDomain::accumulate_iteration_flux() } } -// Compute new estimate of scattering + fission sources in each source region -// based on the flux estimate from the previous iteration. -void FlatSourceDomain::update_neutron_source(double k_eff) +void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } // Add scattering + fission source -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; double scatter_source = 0.0; double fission_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = source_regions_.scalar_flux_old(sr, g_in); + double scalar_flux = srh.scalar_flux_old(g_in); double sigma_s = sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; @@ -154,18 +126,30 @@ void FlatSourceDomain::update_neutron_source(double k_eff) scatter_source += sigma_s * scalar_flux; fission_source += nu_sigma_f * scalar_flux * chi; } - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_source + fission_source * inverse_k_eff) / sigma_t; } } // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } +} + +// Compute new estimate of scattering + fission sources in each source region +// based on the flux estimate from the previous iteration. +void FlatSourceDomain::update_all_neutron_sources() +{ + simulation::time_update_src.start(); + +#pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions(); sr++) { + SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); + update_single_neutron_source(srh); + } simulation::time_update_src.stop(); } @@ -320,7 +304,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() // Generates new estimate of k_eff based on the differences between this // iteration's estimate of the scalar flux and the last iteration's estimate. -double FlatSourceDomain::compute_k_eff(double k_eff_old) const +void FlatSourceDomain::compute_k_eff() { double fission_rate_old = 0; double fission_rate_new = 0; @@ -365,7 +349,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const p[sr] = sr_fission_source_new; } - double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old); double H = 0.0; // defining an inverse sum for better performance @@ -385,7 +369,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Adds entropy value to shared entropy vector in openmc namespace. simulation::entropy.push_back(H); - return k_eff_new; + k_eff_ = k_eff_new; } // This function is responsible for generating a mapping between random @@ -652,7 +636,6 @@ void FlatSourceDomain::random_ray_tally() "random ray mode."); break; } - // Apply score to the appropriate tally bin Tally& tally {*model::tallies[task.tally_idx]}; #pragma omp atomic @@ -726,21 +709,21 @@ void FlatSourceDomain::output_to_vtk() const print_plot(); // Outer loop over plots - for (int p = 0; p < model::plots.size(); p++) { + for (int plt = 0; plt < model::plots.size(); plt++) { // Get handle to OpenMC plot object and extract params - Plot* openmc_plot = dynamic_cast(model::plots[p].get()); + Plot* openmc_plot = dynamic_cast(model::plots[plt].get()); // Random ray plots only support voxel plots if (!openmc_plot) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } else if (openmc_plot->type_ != Plot::PlotType::voxel) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } @@ -794,23 +777,11 @@ void FlatSourceDomain::output_to_vtk() const continue; } - int i_cell = p.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + p.cell_instance(); - if (RandomRay::mesh_subdivision_enabled_) { - int mesh_idx = base_source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(p.r()); - } - SourceRegionKey sr_key {sr, mesh_bin}; - auto it = source_region_map_.find(sr_key); - if (it != source_region_map_.end()) { - sr = it->second; - } else { - sr = -1; - } + SourceRegionKey sr_key = lookup_source_region_key(p); + int64_t sr = -1; + auto it = source_region_map_.find(sr_key); + if (it != source_region_map_.end()) { + sr = it->second; } voxel_indices[z * Ny * Nx + y * Nx + x] = sr; @@ -967,13 +938,17 @@ void FlatSourceDomain::output_to_vtk() const } void FlatSourceDomain::apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh) + int src_idx, SourceRegionHandle& srh) { - srh.external_source_present() = 1; - + auto s = model::external_sources[src_idx].get(); + auto is = dynamic_cast(s); + auto discrete = dynamic_cast(is->energy()); + double strength_factor = is->strength(); const auto& discrete_energies = discrete->x(); const auto& discrete_probs = discrete->prob(); + srh.external_source_present() = 1; + for (int i = 0; i < discrete_energies.size(); i++) { int g = data::mg.get_group_index(discrete_energies[i]); srh.external_source(g) += discrete_probs[i] * strength_factor; @@ -981,8 +956,7 @@ void FlatSourceDomain::apply_external_source_to_source_region( } void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, - const vector& instances) + int src_idx, int target_material_id, const vector& instances) { Cell& cell = *model::cells[i_cell]; @@ -1000,16 +974,13 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, if (target_material_id == C_NONE || cell_material_id == target_material_id) { int64_t source_region = source_region_offsets_[i_cell] + j; - SourceRegionHandle srh = - source_regions_.get_source_region_handle(source_region); - apply_external_source_to_source_region(discrete, strength_factor, srh); + external_volumetric_source_map_[source_region].push_back(src_idx); } } } void FlatSourceDomain::apply_external_source_to_cell_and_children( - int32_t i_cell, Discrete* discrete, double strength_factor, - int32_t target_material_id) + int32_t i_cell, int src_idx, int32_t target_material_id) { Cell& cell = *model::cells[i_cell]; @@ -1017,14 +988,14 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children( vector instances(cell.n_instances()); std::iota(instances.begin(), instances.end(), 0); apply_external_source_to_cell_instances( - i_cell, discrete, strength_factor, target_material_id, instances); + i_cell, src_idx, target_material_id, instances); } else if (target_material_id == C_NONE) { std::unordered_map> cell_instance_list = cell.get_contained_cells(0, nullptr); for (const auto& pair : cell_instance_list) { int32_t i_child_cell = pair.first; - apply_external_source_to_cell_instances(i_child_cell, discrete, - strength_factor, target_material_id, pair.second); + apply_external_source_to_cell_instances( + i_child_cell, src_idx, target_material_id, pair.second); } } } @@ -1070,36 +1041,17 @@ void FlatSourceDomain::convert_external_sources() "point source at {}", sp->r())); } - int i_cell = gs.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); + SourceRegionKey key = lookup_source_region_key(gs); - if (RandomRay::mesh_subdivision_enabled_) { - // If mesh subdivision is enabled, we need to determine which subdivided - // mesh bin the point source coordinate is in as well - int mesh_idx = source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r()); - } - // With the source region and mesh bin known, we can use the - // accompanying SourceRegionKey as a key into a map that stores the - // corresponding external source index for the point source. Notably, we - // do not actually apply the external source to any source regions here, - // as if mesh subdivision is enabled, they haven't actually been - // discovered & initilized yet. When discovered, they will read from the - // point_source_map to determine if there are any point source terms - // that should be applied. - SourceRegionKey key {sr, mesh_bin}; - point_source_map_[key] = es; - } else { - // If we are not using mesh subdivision, we can apply the external - // source directly to the source region as we do for volumetric domain - // constraint sources. - SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); - apply_external_source_to_source_region(energy, strength_factor, srh); - } + // With the source region and mesh bin known, we can use the + // accompanying SourceRegionKey as a key into a map that stores the + // corresponding external source index for the point source. Notably, we + // do not actually apply the external source to any source regions here, + // as if mesh subdivision is enabled, they haven't actually been + // discovered & initilized yet. When discovered, they will read from the + // external_source_map to determine if there are any external source + // terms that should be applied. + external_point_source_map_[key].push_back(es); } else { // If not a point source, then use the volumetric domain constraints to @@ -1107,42 +1059,25 @@ void FlatSourceDomain::convert_external_sources() if (is->domain_type() == Source::DomainType::MATERIAL) { for (int32_t material_id : domain_ids) { for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, material_id); + apply_external_source_to_cell_and_children(i_cell, es, material_id); } } } else if (is->domain_type() == Source::DomainType::CELL) { for (int32_t cell_id : domain_ids) { int32_t i_cell = model::cell_map[cell_id]; - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } else if (is->domain_type() == Source::DomainType::UNIVERSE) { for (int32_t universe_id : domain_ids) { int32_t i_universe = model::universe_map[universe_id]; Universe& universe = *model::universes[i_universe]; for (int32_t i_cell : universe.cells_) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } } } } // End loop over external sources - -// Divide the fixed source term by sigma t (to save time when applying each -// iteration) -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - source_regions_.external_source(sr, g) /= sigma_t; - } - } } void FlatSourceDomain::flux_swap() @@ -1159,13 +1094,23 @@ void FlatSourceDomain::flatten_xs() const int a = 0; n_materials_ = data::mg.macro_xs_.size(); - for (auto& m : data::mg.macro_xs_) { + for (int i = 0; i < n_materials_; i++) { + auto& m = data::mg.macro_xs_[i]; for (int g_out = 0; g_out < negroups_; g_out++) { if (m.exists_in_model) { double sigma_t = m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a); sigma_t_.push_back(sigma_t); + if (sigma_t < MINIMUM_MACRO_XS) { + Material* mat = model::materials[i].get(); + warning(fmt::format( + "Material \"{}\" (id: {}) has a group {} total cross section " + "({:.3e}) below the minimum threshold " + "({:.3e}). Material will be treated as pure void.", + mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS)); + } + double nu_sigma_f = m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a); nu_sigma_f_.push_back(nu_sigma_f); @@ -1206,7 +1151,7 @@ void FlatSourceDomain::flatten_xs() } } -void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) +void FlatSourceDomain::set_adjoint_sources() { // Set the adjoint external source to 1/forward_flux. If the forward flux is // negative, zero, or extremely close to zero, set the adjoint source to zero, @@ -1220,7 +1165,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) double max_flux = 0.0; #pragma omp parallel for reduction(max : max_flux) for (int64_t se = 0; se < n_source_elements(); se++) { - double flux = forward_flux[se]; + double flux = source_regions_.scalar_flux_final(se); if (flux > max_flux) { max_flux = flux; } @@ -1230,7 +1175,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { for (int g = 0; g < negroups_; g++) { - double flux = forward_flux[sr * negroups_ + g]; + double flux = source_regions_.scalar_flux_final(sr, g); if (flux <= ZERO_FLUX_CUTOFF * max_flux) { source_regions_.external_source(sr, g) = 0.0; } else { @@ -1239,6 +1184,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) if (flux > 0.0) { source_regions_.external_source_present(sr) = 1; } + source_regions_.scalar_flux_final(sr, g) = 0.0; } } @@ -1265,7 +1211,6 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) source_regions_.external_source_present(sr) = 0; } } - // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for @@ -1326,13 +1271,14 @@ void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell, if ((target_material_id == C_NONE && !is_target_void) || cell_material_id == target_material_id) { int64_t sr = source_region_offsets_[i_cell] + j; - if (source_regions_.mesh(sr) != C_NONE) { - // print out the source region that is broken: + // Check if the key is already present in the mesh_map_ + if (mesh_map_.find(sr) != mesh_map_.end()) { fatal_error(fmt::format("Source region {} already has mesh idx {} " "applied, but trying to apply mesh idx {}", - sr, source_regions_.mesh(sr), mesh_idx)); + sr, mesh_map_[sr], mesh_idx)); } - source_regions_.mesh(sr) = mesh_idx; + // If the SR has not already been assigned, then we can write to it + mesh_map_[sr] = mesh_idx; } } } @@ -1402,18 +1348,9 @@ void FlatSourceDomain::apply_meshes() } } -void FlatSourceDomain::prepare_base_source_regions() -{ - std::swap(source_regions_, base_source_regions_); - source_regions_.negroups() = base_source_regions_.negroups(); - source_regions_.is_linear() = base_source_regions_.is_linear(); -} - SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u) + SourceRegionKey sr_key, Position r, Direction u) { - SourceRegionKey sr_key {sr, mesh_bin}; - // Case 1: Check if the source region key is already present in the permanent // map. This is the most common condition, as any source region visited in a // previous power iteration will already be present in the permanent map. If @@ -1475,9 +1412,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( gs.r() = r + TINY_BIT * u; gs.u() = {1.0, 0.0, 0.0}; exhaustive_find_cell(gs); - int gs_i_cell = gs.lowest_coord().cell(); - int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance(); - if (sr_found != sr) { + int64_t sr_found = lookup_base_source_region_idx(gs); + if (sr_found != sr_key.base_source_region_id) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1485,9 +1421,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // Sanity check on mesh bin - int mesh_idx = base_source_regions_.mesh(sr); + int mesh_idx = lookup_mesh_idx(sr_key.base_source_region_id); if (mesh_idx == C_NONE) { - if (mesh_bin != 0) { + if (sr_key.mesh_bin != 0) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1496,7 +1432,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } else { Mesh* mesh = model::meshes[mesh_idx].get(); int bin_found = mesh->get_bin(r + TINY_BIT * u); - if (bin_found != mesh_bin) { + if (bin_found != sr_key.mesh_bin) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1508,26 +1444,60 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // condition only occurs the first time the source region is discovered // (typically in the first power iteration). In this case, we need to handle // creation of the new source region and its storage into the parallel map. - // The new source region is created by copying the base source region, so as - // to inherit material, external source, and some flux properties etc. We - // also pass the base source region id to allow the new source region to - // know which base source region it is derived from. - SourceRegion* sr_ptr = discovered_source_regions_.emplace( - sr_key, {base_source_regions_.get_source_region_handle(sr), sr}); - discovered_source_regions_.unlock(sr_key); + // Additionally, we need to determine the source region's material, initialize + // the starting scalar flux guess, and apply any known external sources. + + // Call the basic constructor for the source region and store in the parallel + // map. + bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + SourceRegion* sr_ptr = + discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; - // Check if the new source region contains a point source and apply it if so - auto it2 = point_source_map_.find(sr_key); - if (it2 != point_source_map_.end()) { - int es = it2->second; - auto s = model::external_sources[es].get(); - auto is = dynamic_cast(s); - auto energy = dynamic_cast(is->energy()); - double strength_factor = is->strength(); - apply_external_source_to_source_region(energy, strength_factor, handle); - int material = handle.material(); - if (material != MATERIAL_VOID) { + // Determine the material + int gs_i_cell = gs.lowest_coord().cell(); + Cell& cell = *model::cells[gs_i_cell]; + int material = cell.material(gs.cell_instance()); + + // If material total XS is extremely low, just set it to void to avoid + // problems with 1/Sigma_t + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + if (sigma_t < MINIMUM_MACRO_XS) { + material = MATERIAL_VOID; + break; + } + } + + handle.material() = material; + + // Store the mesh index (if any) assigned to this source region + handle.mesh() = mesh_idx; + + if (settings::run_mode == RunMode::FIXED_SOURCE) { + // Determine if there are any volumetric sources, and apply them. + // Volumetric sources are specifc only to the base SR idx. + auto it_vol = + external_volumetric_source_map_.find(sr_key.base_source_region_id); + if (it_vol != external_volumetric_source_map_.end()) { + const vector& vol_sources = it_vol->second; + for (int src_idx : vol_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // Determine if there are any point sources, and apply them. + // Point sources are specific to the source region key. + auto it_point = external_point_source_map_.find(sr_key); + if (it_point != external_point_source_map_.end()) { + const vector& point_sources = it_point->second; + for (int src_idx : point_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // Divide external source term by sigma_t + if (material != C_NONE) { for (int g = 0; g < negroups_; g++) { double sigma_t = sigma_t_[material * negroups_ + g]; handle.external_source(g) /= sigma_t; @@ -1535,6 +1505,21 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } } + // Compute the combined source term + update_single_neutron_source(handle); + + // Unlock the parallel map. Note: we may be tempted to release + // this lock earlier, and then just use the source region's lock to protect + // the flux/source initialization stages above. However, the rest of the code + // only protects updates to the new flux and volume fields, and assumes that + // the source is constant for the duration of transport. Thus, using just the + // source region's lock by itself would result in other threads potentially + // reading from the source before it is computed, as they won't use the lock + // when only reading from the SR's source. It would be expensive to protect + // those operations, whereas generating the SR is only done once, so we just + // hold the map's bucket lock until the source region is fully initialized. + discovered_source_regions_.unlock(sr_key); + return handle; } @@ -1620,4 +1605,52 @@ void FlatSourceDomain::apply_transport_stabilization() } } +// Determines the base source region index (i.e., a material filled cell +// instance) that corresponds to a particular location in the geometry. Requires +// that the "gs" object passed in has already been initialized and has called +// find_cell etc. +int64_t FlatSourceDomain::lookup_base_source_region_idx( + const GeometryState& gs) const +{ + int i_cell = gs.lowest_coord().cell(); + int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); + return sr; +} + +// Determines the index of the mesh (if any) that has been applied +// to a particular base source region index. +int FlatSourceDomain::lookup_mesh_idx(int64_t sr) const +{ + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } + return mesh_idx; +} + +// Determines the source region key that corresponds to a particular location in +// the geometry. This takes into account both the base source region index as +// well as the mesh bin if a mesh is applied to this source region for +// subdivision. +SourceRegionKey FlatSourceDomain::lookup_source_region_key( + const GeometryState& gs) const +{ + int64_t sr = lookup_base_source_region_idx(gs); + int64_t mesh_bin = lookup_mesh_bin(sr, gs.r()); + return SourceRegionKey {sr, mesh_bin}; +} + +// Determines the mesh bin that corresponds to a particular base source region +// index and position. +int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const +{ + int mesh_idx = lookup_mesh_idx(sr); + int mesh_bin = 0; + if (mesh_idx != C_NONE) { + mesh_bin = model::meshes[mesh_idx]->get_bin(r); + } + return mesh_bin; +} + } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 81412164ec..e1ad68e3d8 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -34,25 +34,18 @@ void LinearSourceDomain::batch_reset() } } -void LinearSourceDomain::update_neutron_source(double k_eff) +void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - MomentMatrix invM = source_regions_.mom_matrix(sr).inverse(); + // Add scattering + fission source + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; + MomentMatrix invM = srh.mom_matrix().inverse(); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; @@ -64,8 +57,8 @@ void LinearSourceDomain::update_neutron_source(double k_eff) for (int g_in = 0; g_in < negroups_; g_in++) { // Handles for the flat and linear components of the flux - double flux_flat = source_regions_.scalar_flux_old(sr, g_in); - MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in); + double flux_flat = srh.scalar_flux_old(g_in); + MomentArray flux_linear = srh.flux_moments_old(g_in); // Handles for cross sections double sigma_s = @@ -81,7 +74,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff) } // Compute the flat source term - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_flat + fission_flat * inverse_k_eff) / sigma_t; // Compute the linear source terms. In the first 10 iterations when the @@ -91,25 +84,21 @@ void LinearSourceDomain::update_neutron_source(double k_eff) // very small/noisy or have poorly developed spatial moments, so we zero // the source gradients (effectively making this a flat source region // temporarily), so as to improve stability. - if (simulation::current_batch > 10 && - source_regions_.source(sr, g_out) >= 0.0) { - source_regions_.source_gradients(sr, g_out) = + if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) { + srh.source_gradients(g_out) = invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); } else { - source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0}; + srh.source_gradients(g_out) = {0.0, 0.0, 0.0}; } } } + // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -// Add external source to flat source term if in fixed source mode -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } - - simulation::time_update_src.stop(); } void LinearSourceDomain::normalize_scalar_flux_and_volumes( diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 27f674c427..89a91449df 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -237,7 +237,6 @@ double RandomRay::distance_inactive_; double RandomRay::distance_active_; unique_ptr RandomRay::ray_source_; RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT}; -bool RandomRay::mesh_subdivision_enabled_ {false}; RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG}; RandomRay::RandomRay() @@ -336,71 +335,60 @@ void RandomRay::event_advance_ray() void RandomRay::attenuate_flux(double distance, bool is_active, double offset) { - // Determine source region index etc. - int i_cell = lowest_coord().cell(); - - // The base source region is the spatial region index - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); + // Lookup base source region index + int64_t sr = domain_->lookup_base_source_region_idx(*this); // Perform ray tracing across mesh - if (mesh_subdivision_enabled_) { - // Determine the mesh index for the base source region, if any - int mesh_idx = domain_->base_source_regions_.mesh(sr); + // Determine the mesh index for the base source region, if any + int mesh_idx = domain_->lookup_mesh_idx(sr); - if (mesh_idx == C_NONE) { - // If there's no mesh being applied to this cell, then - // we just attenuate the flux as normal, and set - // the mesh bin to 0 - attenuate_flux_inner(distance, is_active, sr, 0, r()); - } else { - // If there is a mesh being applied to this cell, then - // we loop over all the bin crossings and attenuate - // separately. - Mesh* mesh = model::meshes[mesh_idx].get(); - - // We adjust the start and end positions of the ray slightly - // to accomodate for floating point precision issues that tend - // to occur at mesh boundaries that overlap with geometry lattice - // boundaries. - Position start = r() + (offset + TINY_BIT) * u(); - Position end = start + (distance - 2.0 * TINY_BIT) * u(); - double reduced_distance = (end - start).norm(); - - // Ray trace through the mesh and record bins and lengths - mesh_bins_.resize(0); - mesh_fractional_lengths_.resize(0); - mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); - - // Loop over all mesh bins and attenuate flux - for (int b = 0; b < mesh_bins_.size(); b++) { - double physical_length = reduced_distance * mesh_fractional_lengths_[b]; - attenuate_flux_inner( - physical_length, is_active, sr, mesh_bins_[b], start); - start += physical_length * u(); - } - } + if (mesh_idx == C_NONE) { + // If there's no mesh being applied to this cell, then + // we just attenuate the flux as normal, and set + // the mesh bin to 0 + attenuate_flux_inner(distance, is_active, sr, 0, r()); } else { - attenuate_flux_inner(distance, is_active, sr, C_NONE, r()); + // If there is a mesh being applied to this cell, then + // we loop over all the bin crossings and attenuate + // separately. + Mesh* mesh = model::meshes[mesh_idx].get(); + + // We adjust the start and end positions of the ray slightly + // to accomodate for floating point precision issues that tend + // to occur at mesh boundaries that overlap with geometry lattice + // boundaries. + Position start = r() + (offset + TINY_BIT) * u(); + Position end = start + (distance - 2.0 * TINY_BIT) * u(); + double reduced_distance = (end - start).norm(); + + // Ray trace through the mesh and record bins and lengths + mesh_bins_.resize(0); + mesh_fractional_lengths_.resize(0); + mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); + + // Loop over all mesh bins and attenuate flux + for (int b = 0; b < mesh_bins_.size(); b++) { + double physical_length = reduced_distance * mesh_fractional_lengths_[b]; + attenuate_flux_inner( + physical_length, is_active, sr, mesh_bins_[b], start); + start += physical_length * u(); + } } } void RandomRay::attenuate_flux_inner( double distance, bool is_active, int64_t sr, int mesh_bin, Position r) { + SourceRegionKey sr_key {sr, mesh_bin}; SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - srh = domain_->get_subdivided_source_region_handle( - sr, mesh_bin, r, distance, u()); - if (srh.is_numerical_fp_artifact_) { - return; - } - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); + srh = domain_->get_subdivided_source_region_handle(sr_key, r, u()); + if (srh.is_numerical_fp_artifact_) { + return; } switch (source_shape_) { case RandomRaySourceShape::FLAT: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_flat_source_void(srh, distance, is_active, r); } else { attenuate_flux_flat_source(srh, distance, is_active, r); @@ -408,7 +396,7 @@ void RandomRay::attenuate_flux_inner( break; case RandomRaySourceShape::LINEAR: case RandomRaySourceShape::LINEAR_XY: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_linear_source_void(srh, distance, is_active, r); } else { attenuate_flux_linear_source(srh, distance, is_active, r); @@ -439,7 +427,7 @@ void RandomRay::attenuate_flux_flat_source( n_event()++; // Get material - int material = this->material(); + int material = srh.material(); // MOC incoming flux attenuation + source contribution/attenuation equation for (int g = 0; g < negroups_; g++) { @@ -490,7 +478,7 @@ void RandomRay::attenuate_flux_flat_source_void( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping @@ -537,7 +525,7 @@ void RandomRay::attenuate_flux_linear_source( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); Position& centroid = srh.centroid(); Position midpoint = r + u() * (distance / 2.0); @@ -810,27 +798,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) cell_born() = lowest_coord().cell(); } + SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); + SourceRegionHandle srh = + domain_->get_subdivided_source_region_handle(sr_key, r(), u()); + // Initialize ray's starting angular flux to starting location's isotropic // source - int i_cell = lowest_coord().cell(); - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); - - SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - int mesh_idx = domain_->base_source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - Mesh* mesh = model::meshes[mesh_idx].get(); - mesh_bin = mesh->get_bin(r()); - } - srh = - domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u()); - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); - } - if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { angular_flux_[g] = srh.source(g); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b84..d475b2593e 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -47,97 +47,82 @@ void openmc_run_random_ray() if (mpi::master) validate_random_ray_inputs(); - // Declare forward flux so that it can be saved for later adjoint simulation - vector forward_flux; - SourceRegionContainer forward_source_regions; - SourceRegionContainer forward_base_source_regions; - std::unordered_map - forward_source_region_map; + // Initialize Random Ray Simulation Object + RandomRaySimulation sim; - { - // Initialize Random Ray Simulation Object - RandomRaySimulation sim; + // Initialize fixed sources, if present + sim.apply_fixed_sources_and_mesh_domains(); - // Initialize fixed sources, if present - sim.apply_fixed_sources_and_mesh_domains(); + // Begin main simulation timer + simulation::time_total.start(); - // Begin main simulation timer - simulation::time_total.start(); + // Execute random ray simulation + sim.simulate(); - // Execute random ray simulation - sim.simulate(); + // End main simulation timer + simulation::time_total.stop(); - // End main simulation timer - simulation::time_total.stop(); - - // Normalize and save the final forward flux - sim.domain()->serialize_final_fluxes(forward_flux); - - double source_normalization_factor = - sim.domain()->compute_fixed_source_normalization_factor() / - (settings::n_batches - settings::n_inactive); + // Normalize and save the final forward flux + double source_normalization_factor = + sim.domain()->compute_fixed_source_normalization_factor() / + (settings::n_batches - settings::n_inactive); #pragma omp parallel for - for (uint64_t i = 0; i < forward_flux.size(); i++) { - forward_flux[i] *= source_normalization_factor; - } - - forward_source_regions = sim.domain()->source_regions_; - forward_source_region_map = sim.domain()->source_region_map_; - forward_base_source_regions = sim.domain()->base_source_regions_; - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - sim.output_simulation_results(); + for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) { + sim.domain()->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; } + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + sim.output_simulation_results(); + ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (adjoint_needed) { - reset_timers(); - - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; - - if (mpi::master) - header("ADJOINT FLUX SOLVE", 3); - - // Initialize OpenMC general data structures - openmc_simulation_init(); - - // Initialize Random Ray Simulation Object - RandomRaySimulation adjoint_sim; - - // Initialize adjoint fixed sources, if present - adjoint_sim.prepare_fixed_sources_adjoint(forward_flux, - forward_source_regions, forward_base_source_regions, - forward_source_region_map); - - // Transpose scattering matrix - adjoint_sim.domain()->transpose_scattering_matrix(); - - // Swap nu_sigma_f and chi - adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_); - - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation - adjoint_sim.simulate(); - - // End main simulation timer - simulation::time_total.stop(); - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - adjoint_sim.output_simulation_results(); + if (!adjoint_needed) { + return; } + + reset_timers(); + + // Configure the domain for adjoint simulation + FlatSourceDomain::adjoint_ = true; + + if (mpi::master) + header("ADJOINT FLUX SOLVE", 3); + + // Initialize OpenMC general data structures + openmc_simulation_init(); + + sim.domain()->k_eff_ = 1.0; + + // Initialize adjoint fixed sources, if present + sim.prepare_fixed_sources_adjoint(); + + // Transpose scattering matrix + sim.domain()->transpose_scattering_matrix(); + + // Swap nu_sigma_f and chi + sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); + + // Begin main simulation timer + simulation::time_total.start(); + + // Execute random ray simulation + sim.simulate(); + + // End main simulation timer + simulation::time_total.stop(); + + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + sim.output_simulation_results(); } // Enforces restrictions on inputs in random ray mode. While there are @@ -348,7 +333,6 @@ void validate_random_ray_inputs() // when generating weight windows with FW-CADIS and an overlaid mesh. /////////////////////////////////////////////////////////////////// if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR && - RandomRay::mesh_subdivision_enabled_ && variance_reduction::weight_windows.size() > 0) { warning( "Linear sources may result in negative fluxes in small source regions " @@ -366,7 +350,6 @@ void openmc_reset_random_ray() FlatSourceDomain::mesh_domain_map_.clear(); RandomRay::ray_source_.reset(); RandomRay::source_shape_ = RandomRaySourceShape::FLAT; - RandomRay::mesh_subdivision_enabled_ = false; RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } @@ -412,20 +395,11 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() } } -void RandomRaySimulation::prepare_fixed_sources_adjoint( - vector& forward_flux, SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, - std::unordered_map& - forward_source_region_map) +void RandomRaySimulation::prepare_fixed_sources_adjoint() { + domain_->source_regions_.adjoint_reset(); if (settings::run_mode == RunMode::FIXED_SOURCE) { - if (RandomRay::mesh_subdivision_enabled_) { - domain_->source_regions_ = forward_source_regions; - domain_->source_region_map_ = forward_source_region_map; - domain_->base_source_regions_ = forward_base_source_regions; - domain_->source_regions_.adjoint_reset(); - } - domain_->set_adjoint_sources(forward_flux); + domain_->set_adjoint_sources(); } } @@ -445,22 +419,18 @@ void RandomRaySimulation::simulate() simulation::total_weight = 1.0; // Update source term (scattering + fission) - domain_->update_neutron_source(k_eff_); + domain_->update_all_neutron_sources(); - // Reset scalar fluxes, iteration volume tallies, and region hit flags to - // zero + // Reset scalar fluxes, iteration volume tallies, and region hit flags + // to zero domain_->batch_reset(); - // At the beginning of the simulation, if mesh subvivision is in use, we + // At the beginning of the simulation, if mesh subdivision is in use, we // need to swap the main source region container into the base container, // as the main source region container will be used to hold the true // subdivided source regions. The base container will therefore only // contain the external source region information, the mesh indices, // material properties, and initial guess values for the flux/source. - if (RandomRay::mesh_subdivision_enabled_ && - simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) { - domain_->prepare_base_source_regions(); - } // Start timer for transport simulation::time_transport.start(); @@ -476,11 +446,9 @@ void RandomRaySimulation::simulate() simulation::time_transport.stop(); - // If using mesh subdivision, add any newly discovered source regions - // to the main source region container. - if (RandomRay::mesh_subdivision_enabled_) { - domain_->finalize_discovered_source_regions(); - } + // Add any newly discovered source regions to the main source region + // container. + domain_->finalize_discovered_source_regions(); // Normalize scalar flux and update volumes domain_->normalize_scalar_flux_and_volumes( @@ -494,10 +462,10 @@ void RandomRaySimulation::simulate() if (settings::run_mode == RunMode::EIGENVALUE) { // Compute random ray k-eff - k_eff_ = domain_->compute_k_eff(k_eff_); + domain_->compute_k_eff(); // Store random ray k-eff into OpenMC's native k-eff variable - global_tally_tracklength = k_eff_; + global_tally_tracklength = domain_->k_eff_; } // Execute all tallying tasks, if this is an active batch @@ -507,12 +475,6 @@ void RandomRaySimulation::simulate() // estimate domain_->accumulate_iteration_flux(); - // Generate mapping between source regions and tallies - if (!domain_->mapped_all_tallies_ && - !RandomRay::mesh_subdivision_enabled_) { - domain_->convert_source_regions_to_tallies(0); - } - // Use above mapping to contribute FSR flux data to appropriate // tallies domain_->random_ray_tally(); @@ -522,7 +484,7 @@ void RandomRaySimulation::simulate() domain_->flux_swap(); // Check for any obvious insabilities/nans/infs - instability_check(n_hits, k_eff_, avg_miss_rate_); + instability_check(n_hits, domain_->k_eff_, avg_miss_rate_); } // End MPI master work // Finalize the current batch @@ -571,7 +533,7 @@ void RandomRaySimulation::instability_check( } if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) { - fatal_error("Instability detected"); + fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff)); } } } diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 1205b995a6..3b06f0ed09 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -48,7 +48,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } scalar_flux_new_.assign(negroups, 0.0); - source_.resize(negroups); + source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); tally_task_.resize(negroups); @@ -60,25 +60,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } } -SourceRegion::SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr) - : SourceRegion(handle.negroups_, handle.is_linear_) -{ - material_ = handle.material(); - mesh_ = handle.mesh(); - parent_sr_ = parent_sr; - for (int g = 0; g < scalar_flux_new_.size(); g++) { - scalar_flux_old_[g] = handle.scalar_flux_old(g); - source_[g] = handle.source(g); - } - - if (settings::run_mode == RunMode::FIXED_SOURCE) { - external_source_present_ = handle.external_source_present(); - for (int g = 0; g < scalar_flux_new_.size(); g++) { - external_source_[g] = handle.external_source(g); - } - } -} - //============================================================================== // SourceRegionContainer implementation //============================================================================== @@ -259,9 +240,12 @@ void SourceRegionContainer::adjoint_reset() MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(), MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); - std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + } else { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0); + } std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0); - std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0); std::fill(source_.begin(), source_.end(), 0.0f); std::fill(external_source_.begin(), external_source_.end(), 0.0f); std::fill(source_gradients_.begin(), source_gradients_.end(), diff --git a/src/settings.cpp b/src/settings.cpp index 03d42bb9da..325256cdc5 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -126,6 +126,7 @@ SolverType solver_type {SolverType::MONTE_CARLO}; std::unordered_set sourcepoint_batch; std::unordered_set statepoint_batch; double source_rejection_fraction {0.05}; +double free_gas_threshold {400.0}; std::unordered_set source_write_surf_id; int64_t ssw_max_particles; int64_t ssw_max_files; @@ -346,7 +347,6 @@ void get_run_parameters(pugi::xml_node node_base) } FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back( type, domain_id); - RandomRay::mesh_subdivision_enabled_ = true; } } } @@ -652,6 +652,10 @@ void read_settings_xml(pugi::xml_node root) std::stod(get_node_value(root, "source_rejection_fraction")); } + if (check_for_node(root, "free_gas_threshold")) { + free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold")); + } + // Survival biasing if (check_for_node(root, "survival_biasing")) { survival_biasing = get_node_value_bool(root, "survival_biasing"); diff --git a/src/simulation.cpp b/src/simulation.cpp index b9986f4737..f55546c4ce 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -674,8 +674,9 @@ void calculate_work() void initialize_data() { // Determine minimum/maximum energy for incident neutron/photon data - data::energy_max = {INFTY, INFTY}; - data::energy_min = {0.0, 0.0}; + data::energy_max = {INFTY, INFTY, INFTY, INFTY}; + data::energy_min = {0.0, 0.0, 0.0, 0.0}; + for (const auto& nuc : data::nuclides) { if (nuc->grid_.size() >= 1) { int neutron = static_cast(ParticleType::neutron); @@ -703,11 +704,21 @@ void initialize_data() // than the current minimum/maximum if (data::ttb_e_grid.size() >= 1) { int photon = static_cast(ParticleType::photon); + int electron = static_cast(ParticleType::electron); + int positron = static_cast(ParticleType::positron); int n_e = data::ttb_e_grid.size(); + + const std::vector charged = {electron, positron}; + for (auto t : charged) { + data::energy_min[t] = std::exp(data::ttb_e_grid(1)); + data::energy_max[t] = std::exp(data::ttb_e_grid(n_e - 1)); + } + data::energy_min[photon] = - std::max(data::energy_min[photon], std::exp(data::ttb_e_grid(1))); - data::energy_max[photon] = std::min( - data::energy_max[photon], std::exp(data::ttb_e_grid(n_e - 1))); + std::max(data::energy_min[photon], data::energy_min[electron]); + + data::energy_max[photon] = + std::min(data::energy_max[photon], data::energy_max[electron]); } } } diff --git a/src/source.cpp b/src/source.cpp index 12323f7bd7..f6aa665ebd 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -290,6 +290,12 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) } else if (temp_str == "photon") { particle_ = ParticleType::photon; settings::photon_transport = true; + } else if (temp_str == "electron") { + particle_ = ParticleType::electron; + settings::photon_transport = true; + } else if (temp_str == "positron") { + particle_ = ParticleType::positron; + settings::photon_transport = true; } else { fatal_error(std::string("Unknown source particle type: ") + temp_str); } diff --git a/tests/regression_tests/electron_heating/__init__.py b/tests/regression_tests/electron_heating/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/electron_heating/inputs_true.dat b/tests/regression_tests/electron_heating/inputs_true.dat new file mode 100644 index 0000000000..ec8e5a8376 --- /dev/null +++ b/tests/regression_tests/electron_heating/inputs_true.dat @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + fixed source + 10000 + 1 + + + 10000000.0 1.0 + + + + 1000.0 + + + + + heating + + + diff --git a/tests/regression_tests/electron_heating/results_true.dat b/tests/regression_tests/electron_heating/results_true.dat new file mode 100644 index 0000000000..4f54ceaa4d --- /dev/null +++ b/tests/regression_tests/electron_heating/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +1.000000E+07 +1.000000E+14 diff --git a/tests/regression_tests/electron_heating/test.py b/tests/regression_tests/electron_heating/test.py new file mode 100644 index 0000000000..e7a58560c4 --- /dev/null +++ b/tests/regression_tests/electron_heating/test.py @@ -0,0 +1,40 @@ +import pytest +import openmc + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def water_model(): + # Define materals and geometry + water = openmc.Material() + water.add_nuclide("H1", 2.0) + water.add_nuclide("O16", 1.0) + water.set_density("g/cc", 1.0) + sphere = openmc.Sphere(r=1.0, boundary_type="reflective") + sph = openmc.Cell(fill=water, region=-sphere) + geometry = openmc.Geometry([sph]) + source = openmc.IndependentSource( + energy=openmc.stats.delta_function(10.0e6), + particle="electron" + ) + + # Define settings + settings = openmc.Settings() + settings.particles = 10000 + settings.batches = 1 + settings.cutoff = {"energy_photon": 1000.0} + settings.run_mode = "fixed source" + settings.source = source + + # Define tallies + tally = openmc.Tally() + tally.scores = ["heating"] + tallies = openmc.Tallies([tally]) + + return openmc.Model(geometry=geometry, settings=settings, tallies=tallies) + + +def test_electron_heating_calc(water_model): + harness = PyAPITestHarness("statepoint.1.h5", water_model) + harness.main() diff --git a/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat index 30e62a8853..0adfc54884 100644 --- a/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat +++ b/tests/regression_tests/random_ray_adjoint_fixed_source/inputs_true.dat @@ -212,8 +212,8 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True - True + true + true naive diff --git a/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat index cd4e92aa1b..073348c41e 100644 --- a/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat +++ b/tests/regression_tests/random_ray_adjoint_k_eff/inputs_true.dat @@ -85,8 +85,8 @@ -1.26 -1.26 -1 1.26 1.26 1 - True - True + true + true diff --git a/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat index 4b8af76aaf..9f1987f3ac 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/cell/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat index 82fe48b614..b4f57dbfa8 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/material/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat index c4fd06f421..ab91f74e50 100644 --- a/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_domain/universe/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true 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 index 4085b0c7a1..220fa7db64 100644 --- 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 @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear 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 index ff085d2fbb..f8c4430852 100644 --- 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 @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear_xy diff --git a/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat index 12c4d74edb..c84e544fcc 100644 --- a/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_mesh/flat/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat index 07adb1ff55..05c4846e6b 100644 --- a/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_mesh/linear/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat index ab077ae8aa..0c870e1006 100644 --- a/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_normalization/False/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - False + false diff --git a/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat index c4fd06f421..ab91f74e50 100644 --- a/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_normalization/True/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat index a124b5e3c4..0c05a71df3 100644 --- a/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_subcritical/flat/inputs_true.dat @@ -115,7 +115,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - False + false flat diff --git a/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat index cc557afb6e..a67495bf16 100644 --- a/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat +++ b/tests/regression_tests/random_ray_fixed_source_subcritical/linear_xy/inputs_true.dat @@ -115,7 +115,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - False + false linear_xy diff --git a/tests/regression_tests/random_ray_halton_samples/inputs_true.dat b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat index 3f058e45c2..36d5f6f227 100644 --- a/tests/regression_tests/random_ray_halton_samples/inputs_true.dat +++ b/tests/regression_tests/random_ray_halton_samples/inputs_true.dat @@ -85,7 +85,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - True + true halton diff --git a/tests/regression_tests/random_ray_k_eff/inputs_true.dat b/tests/regression_tests/random_ray_k_eff/inputs_true.dat index 33c9cac339..545bd1d457 100644 --- a/tests/regression_tests/random_ray_k_eff/inputs_true.dat +++ b/tests/regression_tests/random_ray_k_eff/inputs_true.dat @@ -85,7 +85,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - True + true diff --git a/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat b/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat index f0822d9f96..98badea18d 100644 --- a/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat +++ b/tests/regression_tests/random_ray_k_eff_mesh/inputs_true.dat @@ -85,7 +85,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - True + true diff --git a/tests/regression_tests/random_ray_linear/linear/inputs_true.dat b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat index 0069965572..a43a66e71c 100644 --- a/tests/regression_tests/random_ray_linear/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_linear/linear/inputs_true.dat @@ -85,7 +85,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - True + true linear 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 index 81527c4793..7f76f2fd1c 100644 --- a/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat +++ b/tests/regression_tests/random_ray_linear/linear_xy/inputs_true.dat @@ -85,7 +85,7 @@ -1.26 -1.26 -1 1.26 1.26 1 - True + true linear_xy diff --git a/tests/regression_tests/random_ray_low_density/__init__.py b/tests/regression_tests/random_ray_low_density/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/random_ray_low_density/inputs_true.dat b/tests/regression_tests/random_ray_low_density/inputs_true.dat new file mode 100644 index 0000000000..ab91f74e50 --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/inputs_true.dat @@ -0,0 +1,244 @@ + + + + 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 + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_low_density/results_true.dat b/tests/regression_tests/random_ray_low_density/results_true.dat new file mode 100644 index 0000000000..a4b3ee1bcd --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +5.973607E-01 +7.155477E-02 +tally 2: +3.206216E-02 +2.063375E-04 +tally 3: +2.096415E-03 +8.804963E-07 diff --git a/tests/regression_tests/random_ray_low_density/test.py b/tests/regression_tests/random_ray_low_density/test.py new file mode 100644 index 0000000000..1b4ffb7818 --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/test.py @@ -0,0 +1,60 @@ +import os + +import numpy as np +import openmc +from openmc.examples import random_ray_three_region_cube + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_low_density(): + model = random_ray_three_region_cube() + + # Rebuild the MGXS library to have a material with very + # low macroscopic cross sections + ebins = [1e-5, 20.0e6] + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + void_sigma_a = 4.0e-6 + void_sigma_s = 3.0e-4 + void_mat_data = openmc.XSdata('void', groups) + void_mat_data.order = 0 + void_mat_data.set_total([void_sigma_a + void_sigma_s]) + void_mat_data.set_absorption([void_sigma_a]) + void_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[void_sigma_s]]]), 0, 3)) + + absorber_sigma_a = 0.75 + absorber_sigma_s = 0.25 + absorber_mat_data = openmc.XSdata('absorber', groups) + absorber_mat_data.order = 0 + absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s]) + absorber_mat_data.set_absorption([absorber_sigma_a]) + absorber_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3)) + + multiplier = 0.0000001 + source_sigma_a = void_sigma_a * multiplier + source_sigma_s = void_sigma_s * multiplier + source_mat_data = openmc.XSdata('source', groups) + source_mat_data.order = 0 + source_mat_data.set_total([source_sigma_a + source_sigma_s]) + source_mat_data.set_absorption([source_sigma_a]) + source_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3)) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas( + [source_mat_data, void_mat_data, absorber_mat_data]) + mg_cross_sections_file.export_to_hdf5() + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat index 82ad979da8..088f803bfa 100644 --- a/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat +++ b/tests/regression_tests/random_ray_point_source_locator/inputs_true.dat @@ -211,7 +211,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/random_ray_point_source_locator/results_true.dat b/tests/regression_tests/random_ray_point_source_locator/results_true.dat index 1785dda574..8c6f358dd3 100644 --- a/tests/regression_tests/random_ray_point_source_locator/results_true.dat +++ b/tests/regression_tests/random_ray_point_source_locator/results_true.dat @@ -1,9 +1,9 @@ tally 1: -2.633900E+00 -2.948207E+00 +2.633923E+00 +2.948228E+00 tally 2: -1.440463E-01 -3.294032E-03 +1.440456E-01 +3.293984E-03 tally 3: 9.425207E-03 1.089748E-05 diff --git a/tests/regression_tests/random_ray_void/flat/inputs_true.dat b/tests/regression_tests/random_ray_void/flat/inputs_true.dat index ea8c22f0b5..aa28e7b68b 100644 --- a/tests/regression_tests/random_ray_void/flat/inputs_true.dat +++ b/tests/regression_tests/random_ray_void/flat/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true flat diff --git a/tests/regression_tests/random_ray_void/linear/inputs_true.dat b/tests/regression_tests/random_ray_void/linear/inputs_true.dat index a089604ef4..e4b2f22fa2 100644 --- a/tests/regression_tests/random_ray_void/linear/inputs_true.dat +++ b/tests/regression_tests/random_ray_void/linear/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear diff --git a/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat index 089534d744..8e8a8ed9b8 100644 --- a/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/hybrid/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true hybrid diff --git a/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat index 56507df045..1e25b97da6 100644 --- a/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/naive/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true naive diff --git a/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat index 0331b562cd..78c1626976 100644 --- a/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator/simulation_averaged/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true simulation_averaged diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat index 343833210e..47a8a71824 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/hybrid/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear hybrid diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat index 50be43eb38..80a9ada4d5 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/naive/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear naive diff --git a/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat b/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat index 3d64ab9783..4f032a62a8 100644 --- a/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat +++ b/tests/regression_tests/random_ray_volume_estimator_linear/simulation_averaged/inputs_true.dat @@ -212,7 +212,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true linear simulation_averaged diff --git a/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat index 448b4b145c..5fa6505ddf 100644 --- a/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis/inputs_true.dat @@ -227,7 +227,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true naive diff --git a/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat index d82aa6fae8..ceb89e6e34 100644 --- a/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis_mesh/flat/inputs_true.dat @@ -227,7 +227,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat b/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat index 9e4b21d27b..c7691e950c 100644 --- a/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat +++ b/tests/regression_tests/weightwindows_fw_cadis_mesh/linear/inputs_true.dat @@ -227,7 +227,7 @@ 0.0 0.0 0.0 30.0 30.0 30.0 - True + true diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index cbe00d680c..eace1976ef 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -49,6 +49,7 @@ ENERGIES = np.logspace(log10(1e-5), log10(2e7), 100) ("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)']}, 1e-5), ("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)'], 'nuclides': ['W186', 'H3']}, 1e-2), ]) +@pytest.mark.flaky(reruns=1) def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts, tolerance): # Determine (n.gamma) reaction rate using initial run sp = model.run() diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 6611e4227e..fe618fd2d6 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -59,16 +59,28 @@ def test_export_to_xml(run_in_tmpdir): s.electron_treatment = 'led' s.write_initial_source = True s.weight_window_checkpoints = {'surface': True, 'collision': False} + source_region_mesh = openmc.RegularMesh() + source_region_mesh.dimension = [2, 2, 2] + source_region_mesh.lower_left = [-2, -2, -2] + source_region_mesh.upper_right = [2, 2, 2] + root_universe = openmc.Universe() s.random_ray = { 'distance_inactive': 10.0, 'distance_active': 100.0, 'ray_source': openmc.IndependentSource( space=openmc.stats.Box((-1., -1., -1.), (1., 1., 1.)) - ) + ), + 'source_region_meshes': [(source_region_mesh, [root_universe])], + 'volume_estimator': 'hybrid', + 'source_shape': 'linear', + 'volume_normalized_flux_tallies': True, + 'adjoint': False, + 'sample_method': 'halton' } s.max_particle_events = 100 s.max_secondaries = 1_000_000 s.source_rejection_fraction = 0.01 + s.free_gas_threshold = 800.0 # Make sure exporting XML works s.export_to_xml() @@ -145,5 +157,18 @@ def test_export_to_xml(run_in_tmpdir): assert s.random_ray['distance_active'] == 100.0 assert s.random_ray['ray_source'].space.lower_left == [-1., -1., -1.] assert s.random_ray['ray_source'].space.upper_right == [1., 1., 1.] + assert 'source_region_meshes' in s.random_ray + assert len(s.random_ray['source_region_meshes']) == 1 + mesh_and_domains = s.random_ray['source_region_meshes'][0] + recovered_mesh = mesh_and_domains[0] + assert recovered_mesh.dimension == (2, 2, 2) + assert recovered_mesh.lower_left == [-2., -2., -2.] + assert recovered_mesh.upper_right == [2., 2., 2.] + assert s.random_ray['volume_estimator'] == 'hybrid' + assert s.random_ray['source_shape'] == 'linear' + assert s.random_ray['volume_normalized_flux_tallies'] + assert not s.random_ray['adjoint'] + assert s.random_ray['sample_method'] == 'halton' assert s.max_secondaries == 1_000_000 assert s.source_rejection_fraction == 0.01 + assert s.free_gas_threshold == 800.0 diff --git a/tests/unit_tests/test_statepoint.py b/tests/unit_tests/test_statepoint.py new file mode 100644 index 0000000000..7ffaf7ec2c --- /dev/null +++ b/tests/unit_tests/test_statepoint.py @@ -0,0 +1,65 @@ +import openmc + + +def test_get_tally_filter_type(run_in_tmpdir): + """Test various ways of retrieving tallies from a StatePoint object.""" + + mat = openmc.Material() + mat.add_nuclide("H1", 1.0) + mat.set_density("g/cm3", 10.0) + + sphere = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(fill=mat, region=-sphere) + geometry = openmc.Geometry([cell]) + + settings = openmc.Settings() + settings.particles = 10 + settings.batches = 2 + settings.run_mode = "fixed source" + + reg_mesh = openmc.RegularMesh().from_domain(cell) + tally1 = openmc.Tally(tally_id=1) + mesh_filter = openmc.MeshFilter(reg_mesh) + tally1.filters = [mesh_filter] + tally1.scores = ["flux"] + + tally2 = openmc.Tally(tally_id=2, name="heating tally") + cell_filter = openmc.CellFilter(cell) + tally2.filters = [cell_filter] + tally2.scores = ["heating"] + + tallies = openmc.Tallies([tally1, tally2]) + model = openmc.Model( + geometry=geometry, materials=[mat], settings=settings, tallies=tallies + ) + + sp_filename = model.run() + + sp = openmc.StatePoint(sp_filename) + + tally_found = sp.get_tally(filter_type=openmc.MeshFilter) + assert tally_found.id == 1 + + tally_found = sp.get_tally(filter_type=openmc.CellFilter) + assert tally_found.id == 2 + + tally_found = sp.get_tally(filters=[mesh_filter]) + assert tally_found.id == 1 + + tally_found = sp.get_tally(filters=[cell_filter]) + assert tally_found.id == 2 + + tally_found = sp.get_tally(scores=["heating"]) + assert tally_found.id == 2 + + tally_found = sp.get_tally(name="heating tally") + assert tally_found.id == 2 + + tally_found = sp.get_tally(name=None) + assert tally_found.id == 1 + + tally_found = sp.get_tally(id=1) + assert tally_found.id == 1 + + tally_found = sp.get_tally(id=2) + assert tally_found.id == 2 diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index 386181f34d..abf143f12a 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -16,6 +16,7 @@ def assert_sample_mean(samples, expected_mean): assert np.abs(expected_mean - samples.mean()) < 4*std_dev +@pytest.mark.flaky(reruns=1) def test_discrete(): x = [0.0, 1.0, 10.0] p = [0.3, 0.2, 0.5] @@ -104,6 +105,7 @@ def test_clip_discrete(): d.clip(5) +@pytest.mark.flaky(reruns=1) def test_uniform(): a, b = 10.0, 20.0 d = openmc.stats.Uniform(a, b) @@ -127,6 +129,7 @@ def test_uniform(): assert_sample_mean(samples, exp_mean) +@pytest.mark.flaky(reruns=1) def test_powerlaw(): a, b, n = 10.0, 100.0, 2.0 d = openmc.stats.PowerLaw(a, b, n) @@ -148,6 +151,7 @@ def test_powerlaw(): assert_sample_mean(samples, exp_mean) +@pytest.mark.flaky(reruns=1) def test_maxwell(): theta = 1.2895e6 d = openmc.stats.Maxwell(theta) @@ -171,6 +175,7 @@ def test_maxwell(): assert samples_2.mean() != samples.mean() +@pytest.mark.flaky(reruns=1) def test_watt(): a, b = 0.965e6, 2.29e-6 d = openmc.stats.Watt(a, b) @@ -194,6 +199,7 @@ def test_watt(): assert_sample_mean(samples, exp_mean) +@pytest.mark.flaky(reruns=1) def test_tabular(): # test linear-linear sampling x = np.array([0.0, 5.0, 7.0, 10.0]) @@ -270,6 +276,7 @@ def test_legendre(): d.to_xml_element('distribution') +@pytest.mark.flaky(reruns=1) def test_mixture(): d1 = openmc.stats.Uniform(0, 5) d2 = openmc.stats.Uniform(3, 7) @@ -425,6 +432,7 @@ def test_point(): assert d.xyz == pytest.approx(p) +@pytest.mark.flaky(reruns=1) def test_normal(): mean = 10.0 std_dev = 2.0 @@ -444,6 +452,7 @@ def test_normal(): assert_sample_mean(samples, mean) +@pytest.mark.flaky(reruns=1) def test_muir(): mean = 10.0 mass = 5.0 @@ -463,6 +472,7 @@ def test_muir(): assert_sample_mean(samples, mean) +@pytest.mark.flaky(reruns=1) def test_combine_distributions(): # Combine two discrete (same data as in test_merge_discrete) x1 = [0.0, 1.0, 10.0] diff --git a/tools/ci/gha-script.sh b/tools/ci/gha-script.sh index c1d2921377..b40238ffb5 100755 --- a/tools/ci/gha-script.sh +++ b/tools/ci/gha-script.sh @@ -15,7 +15,7 @@ if [[ $EVENT == 'y' ]]; then fi # Run unit tests and then regression tests -pytest --cov=openmc -v $args \ +pytest -v $args \ tests/test_matplotlib_import.py \ tests/unit_tests \ tests/regression_tests