diff --git a/docs/source/io_formats/mgxs_library.rst b/docs/source/io_formats/mgxs_library.rst index f7f5387a483..8a26311d1e3 100644 --- a/docs/source/io_formats/mgxs_library.rst +++ b/docs/source/io_formats/mgxs_library.rst @@ -133,6 +133,10 @@ Temperature-dependent data, provided for temperature K. This dataset is optional. This is a 1-D vector if `representation` is "isotropic", or a 3-D vector if `representation` is "angle" with dimensions of [polar][azimuthal][groups]. + When this data is not available, an approximation using the + group energy boundaries is used. For more information see + the particle speed subsection in the multigroup-data section + of the theory manual. **//K/scatter_data/** diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index a66abb3ed40..567d24720e9 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -289,6 +289,47 @@ sections. This allows flexibility for the model to use highly anisotropic scattering information in the water while the fuel can be simulated with linear or even isotropic scattering. +Particle Speed +-------------- + +When using a multigroup representation of cross sections, the particle speed has +meaning only in an average sense. The particle speed is important when modeling +dynamic behavior. OpenMC calculates the particle speed using the inverse +velocity multigroup data if it is available. If such data is not available, +OpenMC uses an approximate velocity using the group energy bounds in the +following way: + +.. math:: + + \frac{1}{v_g} = \int_{E_{\text{min}}^g}^{E_{\text{max}}^g} \frac{1}{v(E)} \frac{k}{E} dE + +Where :math:`E_{\text{min}}^g` and :math:`E_{\text{max}}^g` are the group energy +boundaries for group :math:`g`. :math:`v(E)` is the neutron velocity calculated +using relativistic kinematics, :math:`k` is a normalization constant for the +:math:`\frac{1}{E}` spectrum. + +This equation is valid when inside the group boundaries the neutron spectrum +follows a typical :math:`\frac{1}{E}` slowing down spectrum. This assumption is +widely used when generating fine group neutron cross section data libraries from +continuous energy data. + +The solution to this equation is: + +.. math:: + + v_g = \frac{1}{c \log\left(\frac{E_{\text{max}}^g}{E_{\text{min}}^g}\right)} + \left[ 2(\arctan(k_\text{max}) - \arctan(k_\text{min})) + (k_\text{max}-k_\text{min})) \right] + +where :math:`k_{\text{max}}`, :math:`k_{\text{min}}` are defined by a change of +variables: + +.. math:: + + k = \sqrt{1+\frac{2 m_n}{E}} + +where :math:`E` is the particle kinetic energy in [eV] and :math:`m_n` is the +neutron mass in units of [eV]. + .. _logarithmic mapping technique: https://mcnp.lanl.gov/pdf_files/TechReport_2014_LANL_LA-UR-14-24530_Brown.pdf .. _Hwang: https://doi.org/10.13182/NSE87-A16381 diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index da074f825ee..117ac503d49 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -61,6 +61,8 @@ class MgxsInterface { vector energy_bin_avg_; vector rev_energy_bins_; vector> nuc_temps_; // all available temperatures + vector + default_inverse_velocity_; // approximate default inverse-velocity data }; namespace data { diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 2f6e6196bf1..e75f3785a53 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -39,6 +39,8 @@ class Particle : public ParticleData { double speed() const; + double mass() const; + //! create a secondary particle // //! stores the current phase space attributes of the particle in the diff --git a/openmc/mgxs/groups.py b/openmc/mgxs/groups.py index 8910c7d423d..aea7a6d2995 100644 --- a/openmc/mgxs/groups.py +++ b/openmc/mgxs/groups.py @@ -78,6 +78,7 @@ def group_edges(self): @group_edges.setter def group_edges(self, edges): cv.check_type('group edges', edges, Iterable, Real) + cv.check_increasing('group edges', edges) cv.check_greater_than('number of group edges', len(edges), 1) self._group_edges = np.array(edges) diff --git a/src/mgxs.cpp b/src/mgxs.cpp index a0fe4060ddb..1dc090ab969 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -510,6 +510,8 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, break; case MgxsType::INVERSE_VELOCITY: val = xs_t->inverse_velocity(a, gin); + if (!(val > 0)) + val = data::mg.default_inverse_velocity_[gin]; break; case MgxsType::DECAY_RATE: if (dg != nullptr) { diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 08c64413a87..c2253ce174b 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -237,6 +237,18 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) "library file!"); } + // Calculate approximate default inverse velocity data + for (int i = 0; i < energy_bins_.size() - 1; ++i) { + double e_min = std::max(energy_bins_[i + 1], 1e-5); + double e_max = energy_bins_[i]; + double f = 1.0 / (C_LIGHT * std::log(e_max / e_min)); + double k_max = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_max); + double k_min = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_min); + double inv_v = + f * (2.0 * (std::atan(k_max) - std::atan(k_min)) - (k_max - k_min)); + default_inverse_velocity_.push_back(inv_v); + } + // Close MGXS HDF5 file file_close(file_id); } diff --git a/src/particle.cpp b/src/particle.cpp index 0a063559824..3f7691f36fb 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -47,26 +47,33 @@ double Particle::speed() const { if (settings::run_CE) { // Determine mass in eV/c^2 - double mass; - switch (type().pdg_number()) { - case PDG_NEUTRON: - mass = MASS_NEUTRON_EV; - case PDG_ELECTRON: - case PDG_POSITRON: - mass = MASS_ELECTRON_EV; - default: - mass = this->type().mass() * AMU_EV; - } + double mass = this->mass(); // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<E() * (this->E() + 2 * mass)) / (this->E() + mass); } else { - auto& macro_xs = data::mg.macro_xs_[this->material()]; + auto mat = this->material(); + if (mat == MATERIAL_VOID) + return 1.0 / data::mg.default_inverse_velocity_[this->g()]; + auto& macro_xs = data::mg.macro_xs_[mat]; int macro_t = this->mg_xs_cache().t; int macro_a = macro_xs.get_angle_index(this->u()); - return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr, - nullptr, nullptr, macro_t, macro_a); + return 1.0 / macro_xs.get_xs( + MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a); + } +} + +double Particle::mass() const +{ + switch (type().pdg_number()) { + case PDG_NEUTRON: + return MASS_NEUTRON_EV; + case PDG_ELECTRON: + case PDG_POSITRON: + return MASS_ELECTRON_EV; + default: + return this->type().mass() * AMU_EV; } } diff --git a/tests/unit_tests/test_mg_inverse_velocity.py b/tests/unit_tests/test_mg_inverse_velocity.py new file mode 100644 index 00000000000..b520e48d193 --- /dev/null +++ b/tests/unit_tests/test_mg_inverse_velocity.py @@ -0,0 +1,59 @@ +import openmc +import numpy as np +import pytest + +@pytest.fixture +def one_group_lib(): + groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6]) + xsdata = openmc.XSdata('slab_mat', groups) + xsdata.order = 0 + xsdata.set_total([0.0]) + xsdata.set_absorption([0.0]) + xsdata.set_scatter_matrix([[[0.0]]]) + + mg_library = openmc.MGXSLibrary(groups) + mg_library.add_xsdata(xsdata) + name = 'mgxs.h5' + mg_library.export_to_hdf5(name) + yield name + +@pytest.fixture +def slab_model(one_group_lib): + model = openmc.Model() + mat = openmc.Material(name='slab_material') + mat.set_density('macro', 1.0) + mat.add_macroscopic('slab_mat') + + model.materials = openmc.Materials([mat]) + model.materials.cross_sections = one_group_lib + + x_min = openmc.XPlane(x0=0.0, boundary_type='vacuum') + x_max = openmc.XPlane(x0=10.0, boundary_type='vacuum') + + y_min = openmc.YPlane(y0=-10.0, boundary_type='vacuum') + y_max = openmc.YPlane(y0=10.0, boundary_type='vacuum') + z_min = openmc.ZPlane(z0=-10.0, boundary_type='vacuum') + z_max = openmc.ZPlane(z0=19.0, boundary_type='vacuum') + + cell = openmc.Cell(fill=mat, region=+z_min & -x_max & +y_min & -y_max & +z_min & -z_max) + model.geometry = openmc.Geometry([cell]) + + model.settings = openmc.Settings() + model.settings.energy_mode = 'multi-group' + model.settings.run_mode = 'fixed source' + model.settings.batches = 3 + model.settings.particles = 10 + + source = openmc.IndependentSource() + source.space = openmc.stats.Point((5.0, 0.0, 0.0)) + model.settings.source = source + return model + +def test_inverse_velocity(run_in_tmpdir, slab_model): + tally = openmc.Tally() + tally.scores = ['flux','inverse-velocity'] + slab_model.tallies = [tally] + slab_model.run(apply_tally_results=True) + inverse_velocity = tally.mean.squeeze()[1]/tally.mean.squeeze()[0] + + assert inverse_velocity == pytest.approx(1.6144e-5, rel=1e-4)