diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6ec067a12..a371a3c39 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -96,6 +96,17 @@ When the DAGMC mode is enabled, the OpenMC geometry will be read from the file ``dagmc.h5m``. If a :ref:`geometry.xml ` file is present with ``dagmc`` set to ``true``, it will be ignored. +---------------------------- +```` +---------------------------- + +Determines whether to scale the fission photon yield to account for delayed +photon energy. The photon yields are scaled as (EGP + EGD)/EGP where EGP and EGD +are the prompt and delayed photon components of energy release, respectively, +from MF=1, MT=458 on an ENDF evaluation. + + *Default*: true + -------------------------------- ```` Element -------------------------------- diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 7439d2e34..294197b94 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -71,6 +71,10 @@ Coupling and Multi-physics Coupling of OpenMC and Nek5000 within the MOOSE Framework," *Proc. PHYSOR*, Cancun, Mexico, Apr. 22-26 (2018). +- Sterling Harper, Kord Smith, and Benoit Forget, "Faster Monte Carlo + multiphysics using temperature derivatives," *Proc. PHYSOR*, Cancun, Mexico, + Apr. 22-26 (2018). + - Jun Chen, Liangzhi Cao, Chuanqi Zhao, and Zhouyu Liu, "`Development of Subchannel Code SUBSC for high-fidelity multi-physics coupling application `_", *Energy Procedia*, **127**, diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index a1571ea58..dd9147c22 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -97,7 +97,7 @@ public: std::vector index_inelastic_scatter_; private: - void create_derived(); + void create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons); static int XS_TOTAL; static int XS_ABSORPTION; diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 6570a25f7..e8c1d218c 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -29,6 +29,7 @@ extern bool confidence_intervals; //!< use confidence intervals for results? extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)? extern "C" bool cmfd_run; //!< is a CMFD run? extern "C" bool dagmc; //!< indicator of DAGMC geometry +extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed extern "C" bool entropy_on; //!< calculate Shannon entropy? extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? extern bool material_cell_offsets; //!< create material cells offsets? diff --git a/openmc/settings.py b/openmc/settings.py index f25fd681d..27ca52a38 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -40,6 +40,10 @@ class Settings(object): below which particle type will be killed. dagmc : bool Indicate that a CAD-based DAGMC geometry will be used. + delayed_photon_scaling : bool + Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP + where EGP is the energy release of prompt photons and EGD is the energy + release of delayed photons. electron_treatment : {'led', 'ttb'} Whether to deposit all energy from electrons locally ('led') or create secondary bremsstrahlung photons ('ttb'). @@ -219,6 +223,7 @@ class Settings(object): VolumeCalculation, 'volume calculations') self._create_fission_neutrons = None + self._delayed_photon_scaling = None self._material_cell_offsets = None self._log_grid_bins = None @@ -356,6 +361,10 @@ class Settings(object): def create_fission_neutrons(self): return self._create_fission_neutrons + @property + def delayed_photon_scaling(self): + return self._delayed_photon_scaling + @property def material_cell_offsets(self): return self._material_cell_offsets @@ -691,6 +700,11 @@ class Settings(object): create_fission_neutrons, bool) self._create_fission_neutrons = create_fission_neutrons + @delayed_photon_scaling.setter + def delayed_photon_scaling(self, value): + cv.check_type('delayed photon scaling', value, bool) + self._delayed_photon_scaling = value + @material_cell_offsets.setter def material_cell_offsets(self, value): cv.check_type('material cell offsets', value, bool) @@ -930,6 +944,11 @@ class Settings(object): elem = ET.SubElement(root, "create_fission_neutrons") elem.text = str(self._create_fission_neutrons).lower() + def _create_delayed_photon_scaling_subelement(self, root): + if self._delayed_photon_scaling is not None: + elem = ET.SubElement(root, "delayed_photon_scaling") + elem.text = str(self._delayed_photon_scaling).lower() + def _create_material_cell_offsets_subelement(self, root): if self._material_cell_offsets is not None: elem = ET.SubElement(root, "material_cell_offsets") @@ -1166,6 +1185,11 @@ class Settings(object): if text is not None: self.create_fission_neutrons = text in ('true', '1') + def _delayed_photon_scaling_from_xml_element(self, root): + text = get_text(root, 'delayed_photon_scaling') + if text is not None: + self.delayed_photon_scaling = text in ('true', '1') + def _material_cell_offsets_from_xml_element(self, root): text = get_text(root, 'material_cell_offsets') if text is not None: @@ -1225,6 +1249,7 @@ class Settings(object): self._create_resonance_scattering_subelement(root_element) self._create_volume_calcs_subelement(root_element) self._create_create_fission_neutrons_subelement(root_element) + self._create_delayed_photon_scaling_subelement(root_element) self._create_material_cell_offsets_subelement(root_element) self._create_log_grid_bins_subelement(root_element) self._create_dagmc_subelement(root_element) @@ -1291,6 +1316,7 @@ class Settings(object): settings._ufs_mesh_from_xml_element(root) settings._resonance_scattering_from_xml_element(root) settings._create_fission_neutrons_from_xml_element(root) + settings._delayed_photon_scaling_from_xml_element(root) settings._material_cell_offsets_from_xml_element(root) settings._log_grid_bins_from_xml_element(root) settings._dagmc_from_xml_element(root) diff --git a/src/error.cpp b/src/error.cpp index 83e7e8781..d44e8766e 100644 --- a/src/error.cpp +++ b/src/error.cpp @@ -57,7 +57,7 @@ void output(const std::string& message, std::ostream& out, int indent) while (i_start < length) { if (length - i_start < line_len) { // Remainder of message will fit on line - out << message.substr(i_start) << '\n'; + out << message.substr(i_start) << std::endl; break; } else { diff --git a/src/finalize.cpp b/src/finalize.cpp index d359f3885..5bda587ac 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -69,6 +69,7 @@ int openmc_finalize() settings::check_overlaps = false; settings::confidence_intervals = false; settings::create_fission_neutrons = true; + settings::delayed_photon_scaling = true; settings::electron_treatment = ELECTRON_LED; settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; settings::entropy_on = false; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 86f2ee140..824eabb67 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -259,17 +259,24 @@ Nuclide::Nuclide(hid_t group, const std::vector& temperature, int i_nucl } // Read fission energy release data if present + std::unique_ptr prompt_photons; + std::unique_ptr delayed_photons; if (object_exists(group, "fission_energy_release")) { hid_t fer_group = open_group(group, "fission_energy_release"); fission_q_prompt_ = read_function(fer_group, "q_prompt"); fission_q_recov_ = read_function(fer_group, "q_recoverable"); + + // We need prompt/delayed photon energy release for scaling fission photon + // production + prompt_photons = read_function(fer_group, "prompt_photons"); + delayed_photons = read_function(fer_group, "delayed_photons"); close_group(fer_group); } - this->create_derived(); + this->create_derived(prompt_photons.get(), delayed_photons.get()); } -void Nuclide::create_derived() +void Nuclide::create_derived(const Function1D* prompt_photons, const Function1D* delayed_photons) { for (const auto& grid : grid_) { // Allocate and initialize cross section @@ -294,7 +301,21 @@ void Nuclide::create_derived() auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD); for (int k = 0; k < n; ++k) { double E = grid_[t].energy[k+j]; - pprod[k] += xs[k] * (*p.yield_)(E); + + // For fission, artificially increase the photon yield to account + // for delayed photons + double f = 1.0; + if (settings::delayed_photon_scaling) { + if (is_fission(rx->mt_)) { + if (prompt_photons && delayed_photons) { + double energy_prompt = (*prompt_photons)(E); + double energy_delayed = (*delayed_photons)(E); + f = (energy_prompt + energy_delayed)/(energy_prompt); + } + } + } + + pprod[k] += f * xs[k] * (*p.yield_)(E); } } } diff --git a/src/output.cpp b/src/output.cpp index c8b91c4dc..79bd51253 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -92,7 +92,7 @@ void title() // Write number of OpenMP threads std::cout << " OpenMP Threads | " << omp_get_max_threads() << '\n'; #endif - std::cout << '\n'; + std::cout << std::endl; } //============================================================================== @@ -126,7 +126,7 @@ header(const char* msg, int level) { // Print header based on verbosity level. if (settings::verbosity >= level) - std::cout << '\n' << out << "\n\n"; + std::cout << '\n' << out << "\n" << std::endl; } //============================================================================== diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index c69624c53..3bc95f352 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -14,6 +14,8 @@ element settings { (element energy_positron { xsd:double } | attribute energy_positron { xsd:double })? }? & + element delayed_photon_scaling { xsd:boolean }? & + element electron_treatment { ( "led" | "ttb" ) }? & element energy_grid { ( "nuclide" | "log" | "logarithm" | "logarithmic" | "material-union" | "union" ) }? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 46845b094..4e795419d 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -82,6 +82,11 @@ + + + + + diff --git a/src/settings.cpp b/src/settings.cpp index 59d2479be..a4381d0f3 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -44,6 +44,7 @@ bool cmfd_run {false}; bool confidence_intervals {false}; bool create_fission_neutrons {true}; bool dagmc {false}; +bool delayed_photon_scaling {true}; bool entropy_on {false}; bool legendre_to_tabular {true}; bool material_cell_offsets {true}; @@ -778,6 +779,11 @@ void read_settings_xml() } } + // Check whether to scale fission photon yields + if (check_for_node(root, "delayed_photon_scaling")) { + delayed_photon_scaling = get_node_value_bool(root, "delayed_photon_scaling"); + } + // Check whether material cell offsets should be generated if (check_for_node(root, "material_cell_offsets")) { material_cell_offsets = get_node_value_bool(root, "material_cell_offsets"); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 51eae2620..12a2bb123 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1267,6 +1267,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, break; default: + // The default block is really only meant for redundant neutron reactions + // (e.g. 444, 901) + if (p->type_ != Particle::Type::neutron) continue; + if (tally.estimator_ == ESTIMATOR_ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need // to check if it matches the MT number of the event diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 1821bea92..31eab718b 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -1,16 +1,16 @@ - - - - - - - + + + + + + + - + @@ -38,7 +38,7 @@ - 9 + 1 photon diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 1f645ebb9..4ed0a7ba7 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -1,90 +1,72 @@ -from math import pi - -import numpy as np import openmc +import pytest from tests.testing_harness import PyAPITestHarness -class SourceTestHarness(PyAPITestHarness): - def _build_inputs(self): - mat = openmc.Material() - mat.set_density('g/cm3', 2.6989) - mat.add_nuclide('Al27', 1.0) - materials = openmc.Materials([mat]) - materials.export_to_xml() +@pytest.fixture +def model(): + model = openmc.model.Model() + mat = openmc.Material() + mat.set_density('g/cm3', 2.6989) + mat.add_nuclide('Al27', 1.0) + model.materials.append(mat) - cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum') - x_plane_left = openmc.XPlane(-1.0, 'vacuum') - x_plane_center = openmc.XPlane(1.0) - x_plane_right = openmc.XPlane(1.0e9, 'vacuum') + cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum') + x_plane_left = openmc.XPlane(-1.0, 'vacuum') + x_plane_center = openmc.XPlane(1.0) + x_plane_right = openmc.XPlane(1.0e9, 'vacuum') - inner_cyl_left = openmc.Cell() - inner_cyl_right = openmc.Cell() - outer_cyl = openmc.Cell() + inner_cyl_left = openmc.Cell() + inner_cyl_right = openmc.Cell() + outer_cyl = openmc.Cell() - inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center - inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right - outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right) - inner_cyl_right.fill = mat - geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl]) - geometry.export_to_xml() + inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center + inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right + outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right) + inner_cyl_right.fill = mat + model.geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl]) - source = openmc.Source() - source.space = openmc.stats.Point((0, 0, 0)) - source.angle = openmc.stats.Monodirectional() - source.energy = openmc.stats.Discrete([14.0e6], [1.0]) - source.particle = 'neutron' + source = openmc.Source() + source.space = openmc.stats.Point((0, 0, 0)) + source.angle = openmc.stats.Monodirectional() + source.energy = openmc.stats.Discrete([14.0e6], [1.0]) + source.particle = 'neutron' - settings = openmc.Settings() - settings.particles = 10000 - settings.run_mode = 'fixed source' - settings.batches = 1 - settings.photon_transport = True - settings.electron_treatment = 'ttb' - settings.cutoff = {'energy_photon' : 1000.0} - settings.source = source - settings.export_to_xml() + model.settings.particles = 10000 + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.photon_transport = True + model.settings.electron_treatment = 'ttb' + model.settings.cutoff = {'energy_photon' : 1000.0} + model.settings.source = source - surface_filter = openmc.SurfaceFilter(cyl) - particle_filter = openmc.ParticleFilter('photon') - current_tally = openmc.Tally() - current_tally.filters = [surface_filter, particle_filter] - current_tally.scores = ['current'] - tally_tracklength = openmc.Tally() - tally_tracklength.filters = [particle_filter] - tally_tracklength.scores = ['total', 'heating'] - tally_tracklength.nuclides = ['Al27', 'total'] - tally_tracklength.estimator = 'tracklength' - tally_collision = openmc.Tally() - tally_collision.filters = [particle_filter] - tally_collision.scores = ['total', 'heating'] - tally_collision.nuclides = ['Al27', 'total'] - tally_collision.estimator = 'collision' - tally_analog = openmc.Tally() - tally_analog.filters = [particle_filter] - tally_analog.scores = ['total', 'heating'] - tally_analog.nuclides = ['Al27', 'total'] - tally_analog.estimator = 'analog' - tallies = openmc.Tallies([current_tally, tally_tracklength, - tally_collision, tally_analog]) - tallies.export_to_xml() + surface_filter = openmc.SurfaceFilter(cyl) + particle_filter = openmc.ParticleFilter('photon') + current_tally = openmc.Tally() + current_tally.filters = [surface_filter, particle_filter] + current_tally.scores = ['current'] + tally_tracklength = openmc.Tally() + tally_tracklength.filters = [particle_filter] + tally_tracklength.scores = ['total', 'heating'] + tally_tracklength.nuclides = ['Al27', 'total'] + tally_tracklength.estimator = 'tracklength' + tally_collision = openmc.Tally() + tally_collision.filters = [particle_filter] + tally_collision.scores = ['total', 'heating'] + tally_collision.nuclides = ['Al27', 'total'] + tally_collision.estimator = 'collision' + tally_analog = openmc.Tally() + tally_analog.filters = [particle_filter] + tally_analog.scores = ['total', 'heating'] + tally_analog.nuclides = ['Al27', 'total'] + tally_analog.estimator = 'analog' + model.tallies.extend([current_tally, tally_tracklength, + tally_collision, tally_analog]) - def _get_results(self): - with openmc.StatePoint(self._sp_name) as sp: - outstr = '' - for i, tally_ind in enumerate(sp.tallies): - tally = sp.tallies[tally_ind] - results = np.zeros((tally.sum.size * 2, )) - results[0::2] = tally.sum.ravel() - results[1::2] = tally.sum_sq.ravel() - results = ['{0:12.6E}'.format(x) for x in results] - - outstr += 'tally {}:\n'.format(i + 1) - outstr += '\n'.join(results) + '\n' - return outstr + return model -def test_photon_production(): - harness = SourceTestHarness('statepoint.1.h5') +def test_photon_production(model): + harness = PyAPITestHarness('statepoint.1.h5', model) harness.main() diff --git a/tests/regression_tests/photon_production_fission/__init__.py b/tests/regression_tests/photon_production_fission/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/photon_production_fission/inputs_true.dat b/tests/regression_tests/photon_production_fission/inputs_true.dat new file mode 100644 index 000000000..1c496b2f2 --- /dev/null +++ b/tests/regression_tests/photon_production_fission/inputs_true.dat @@ -0,0 +1,49 @@ + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 2 + + + 0 0 0 + + + true + + + + + neutron photon + + + 1 + U235 total + fission heating heating-local + tracklength + + + 1 + U235 total + fission heating heating-local + collision + + + 1 + U235 total + fission heating heating-local + analog + + diff --git a/tests/regression_tests/photon_production_fission/results_true.dat b/tests/regression_tests/photon_production_fission/results_true.dat new file mode 100644 index 000000000..ca4517b4f --- /dev/null +++ b/tests/regression_tests/photon_production_fission/results_true.dat @@ -0,0 +1,77 @@ +k-combined: +2.223956E+00 2.598313E-03 +tally 1: +2.600054E+00 +2.253433E+00 +4.390171E+08 +6.424554E+16 +0.000000E+00 +0.000000E+00 +2.600054E+00 +2.253433E+00 +4.390171E+08 +6.424554E+16 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.764667E+07 +7.568912E+14 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.764667E+07 +7.568912E+14 +0.000000E+00 +0.000000E+00 +tally 2: +2.602512E+00 +2.258176E+00 +4.394176E+08 +6.437709E+16 +0.000000E+00 +0.000000E+00 +2.602512E+00 +2.258176E+00 +4.394176E+08 +6.437709E+16 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.749088E+07 +7.519602E+14 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.749088E+07 +7.519602E+14 +0.000000E+00 +0.000000E+00 +tally 3: +2.663535E+00 +2.364845E+00 +4.394176E+08 +6.437709E+16 +0.000000E+00 +0.000000E+00 +2.663535E+00 +2.364845E+00 +4.394176E+08 +6.437709E+16 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.729796E+07 +7.458273E+14 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.743827E+07 +7.502589E+14 +0.000000E+00 +0.000000E+00 diff --git a/tests/regression_tests/photon_production_fission/test.py b/tests/regression_tests/photon_production_fission/test.py new file mode 100644 index 000000000..f013c9f80 --- /dev/null +++ b/tests/regression_tests/photon_production_fission/test.py @@ -0,0 +1,48 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.model.Model() + mat = openmc.Material() + mat.set_density('g/cm3', 10.0) + mat.add_nuclide('U235', 1.0) + model.materials.append(mat) + + sph = openmc.Sphere(r=100.0, boundary_type='reflective') + cell = openmc.Cell(fill=mat, region=-sph) + model.geometry = openmc.Geometry([cell]) + + model.settings.particles = 1000 + model.settings.batches = 5 + model.settings.inactive = 2 + model.settings.photon_transport = True + model.settings.source = openmc.Source(space=openmc.stats.Point((0, 0, 0))) + + particle_filter = openmc.ParticleFilter(['neutron', 'photon']) + tally_tracklength = openmc.Tally() + tally_tracklength.filters = [particle_filter] + tally_tracklength.scores = ['fission', 'heating', 'heating-local'] + tally_tracklength.nuclides = ['U235', 'total'] + tally_tracklength.estimator = 'tracklength' + tally_collision = openmc.Tally() + tally_collision.filters = [particle_filter] + tally_collision.scores = ['fission', 'heating', 'heating-local'] + tally_collision.nuclides = ['U235', 'total'] + tally_collision.estimator = 'collision' + tally_analog = openmc.Tally() + tally_analog.filters = [particle_filter] + tally_analog.scores = ['fission', 'heating', 'heating-local'] + tally_analog.nuclides = ['U235', 'total'] + tally_analog.estimator = 'analog' + model.tallies.extend([tally_tracklength, tally_collision, tally_analog]) + + return model + + +def test_photon_production(model): + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main()