From c39990accb6d0377fc05f004b0809d08e7a7f384 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Apr 2017 10:59:31 -0500 Subject: [PATCH 1/3] Fix bug tallying individual reaction rates (e.g. n,gamma) with multipole on --- openmc/settings.py | 5 +- src/tally.F90 | 56 ++++++--- tests/test_multipole/inputs_true.dat | 20 ++-- tests/test_multipole/results_true.dat | 31 +++++ tests/test_multipole/test_multipole.py | 150 +++++++++++-------------- 5 files changed, 145 insertions(+), 117 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index bb04defcd8..b1ff4c0c65 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -997,7 +997,10 @@ class Settings(object): for key, value in sorted(self.temperature.items()): element = ET.SubElement(root, "temperature_{}".format(key)) - element.text = str(value) + if isinstance(value, bool): + element.text = str(value).lower() + else: + element.text = str(value) def _create_threads_subelement(self, root): if self._threads is not None: diff --git a/src/tally.F90 b/src/tally.F90 index ebca2bb5dc..96792f76c6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1169,16 +1169,28 @@ contains ! Retrieve temperature and energy grid index and interpolation ! factor i_temp = micro_xs(i_nuclide) % index_temp - i_energy = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor + if (i_temp > 0) then + i_energy = micro_xs(i_nuclide) % index_grid + f = micro_xs(i_nuclide) % interp_factor - associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) - if (i_energy >= xs % threshold) then - score = ((ONE - f) * xs % value(i_energy - & - xs % threshold + 1) + f * xs % value(i_energy - & - xs % threshold + 2)) * atom_density * flux + associate (xs => nuclides(i_nuclide) % reactions(m) % xs(i_temp)) + if (i_energy >= xs % threshold) then + score = ((ONE - f) * xs % value(i_energy - & + xs % threshold + 1) + f * xs % value(i_energy - & + xs % threshold + 2)) * atom_density * flux + end if + end associate + else + ! This block is reached if multipole is turned on and we're in + ! the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + if (score_bin == N_GAMMA) then + score = (micro_xs(i_nuclide) % absorption - & + micro_xs(i_nuclide) % fission) * atom_density * flux + else + score = ZERO end if - end associate + end if end if else @@ -1195,16 +1207,28 @@ contains ! Retrieve temperature and energy grid index and interpolation ! factor i_temp = micro_xs(i_nuc) % index_temp - i_energy = micro_xs(i_nuc) % index_grid - f = micro_xs(i_nuc) % interp_factor + if (i_temp > 0) then + i_energy = micro_xs(i_nuc) % index_grid + f = micro_xs(i_nuc) % interp_factor - associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp)) - if (i_energy >= xs % threshold) then - score = score + ((ONE - f) * xs % value(i_energy - & - xs % threshold + 1) + f * xs % value(i_energy - & - xs % threshold + 2)) * atom_density_ * flux + associate (xs => nuclides(i_nuc) % reactions(m) % xs(i_temp)) + if (i_energy >= xs % threshold) then + score = score + ((ONE - f) * xs % value(i_energy - & + xs % threshold + 1) + f * xs % value(i_energy - & + xs % threshold + 2)) * atom_density_ * flux + end if + end associate + else + ! This block is reached if multipole is turned on and we're + ! in the resolved range. For (n,gamma), use absorption - + ! fission. For everything else, assume it's zero. + if (score_bin == N_GAMMA) then + score = (micro_xs(i_nuc) % absorption - micro_xs(i_nuc) & + % fission) * atom_density_ * flux + else + score = ZERO end if - end associate + end if end if end do end if diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 4163a00e3c..fdedc5425f 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -43,19 +43,13 @@ -1 -1 -1 1 1 1 - True + true 1000 - - - 0 0 0 - 7 7 - 400 400 - - - 0 0 0 - 7 7 - 400 400 - - + + + U235 O16 total + total fission (n,gamma) elastic (n,p) + + diff --git a/tests/test_multipole/results_true.dat b/tests/test_multipole/results_true.dat index 5418a212af..c2111bff47 100644 --- a/tests/test_multipole/results_true.dat +++ b/tests/test_multipole/results_true.dat @@ -1,5 +1,36 @@ k-combined: 1.363786E+00 1.103929E-02 +tally 1: +3.960375E+00 +3.138356E+00 +2.875799E+00 +1.655329E+00 +5.521904E-01 +6.105083E-02 +4.617974E-01 +4.274130E-02 +0.000000E+00 +0.000000E+00 +2.254165E+01 +1.016444E+02 +0.000000E+00 +0.000000E+00 +6.839351E-04 +9.359599E-08 +2.251829E+01 +1.014341E+02 +4.903575E-05 +1.199780E-09 +3.595002E+02 +2.586079E+04 +2.875799E+00 +1.655329E+00 +2.174747E+00 +9.465589E-01 +3.543564E+02 +2.512619E+04 +4.903575E-05 +1.199780E-09 Cell ID = 11 Name = diff --git a/tests/test_multipole/test_multipole.py b/tests/test_multipole/test_multipole.py index 1ab22f6e69..294185312a 100644 --- a/tests/test_multipole/test_multipole.py +++ b/tests/test_multipole/test_multipole.py @@ -4,94 +4,69 @@ import sys sys.path.insert(0, os.pardir) from testing_harness import TestHarness, PyAPITestHarness import openmc +import openmc.model + + +def make_model(): + model = openmc.model.Model() + + # Materials + moderator = openmc.Material(material_id=1) + moderator.set_density('g/cc', 1.0) + moderator.add_nuclide('H1', 2.0) + moderator.add_nuclide('O16', 1.0) + moderator.add_s_alpha_beta('c_H_in_H2O') + + dense_fuel = openmc.Material(material_id=2) + dense_fuel.set_density('g/cc', 4.5) + dense_fuel.add_nuclide('U235', 1.0) + + model.materials += [moderator, dense_fuel] + + # Geometry + c1 = openmc.Cell(cell_id=1, fill=moderator) + mod_univ = openmc.Universe(universe_id=1, cells=(c1,)) + + r0 = openmc.ZCylinder(R=0.3) + c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0) + c11.temperature = [500, 0, 700, 800] + c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0) + fuel_univ = openmc.Universe(universe_id=11, cells=(c11, c12)) + + lat = openmc.RectLattice(lattice_id=101) + lat.dimension = [2, 2] + lat.lower_left = [-2.0, -2.0] + lat.pitch = [2.0, 2.0] + lat.universes = [[fuel_univ]*2]*2 + lat.outer = mod_univ + + x0 = openmc.XPlane(x0=-3.0) + x1 = openmc.XPlane(x0=3.0) + y0 = openmc.YPlane(y0=-3.0) + y1 = openmc.YPlane(y0=3.0) + for s in [x0, x1, y0, y1]: + s.boundary_type = 'reflective' + c101 = openmc.Cell(cell_id=101, fill=lat, region=+x0 & -x1 & +y0 & -y1) + model.geometry.root_universe = openmc.Universe(universe_id=0, cells=(c101,)) + + # Settings + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + model.settings.source = openmc.Source(space=openmc.stats.Box( + [-1, -1, -1], [1, 1, 1])) + model.settings.temperature = {'tolerance': 1000, 'multipole': True} + + # Tallies + tally = openmc.Tally() + tally.nuclides = ['U235', 'O16', 'total'] + tally.scores = ['total', 'fission', '(n,gamma)', 'elastic', '(n,p)'] + model.tallies.append(tally) + + return model + class MultipoleTestHarness(PyAPITestHarness): - def _build_inputs(self): - #################### - # Materials - #################### - - moderator = openmc.Material(material_id=1) - moderator.set_density('g/cc', 1.0) - moderator.add_nuclide('H1', 2.0) - moderator.add_nuclide('O16', 1.0) - moderator.add_s_alpha_beta('c_H_in_H2O') - - dense_fuel = openmc.Material(material_id=2) - dense_fuel.set_density('g/cc', 4.5) - dense_fuel.add_nuclide('U235', 1.0) - - mats_file = openmc.Materials([moderator, dense_fuel]) - mats_file.export_to_xml() - - #################### - # Geometry - #################### - - c1 = openmc.Cell(cell_id=1, fill=moderator) - mod_univ = openmc.Universe(universe_id=1, cells=(c1,)) - - r0 = openmc.ZCylinder(R=0.3) - c11 = openmc.Cell(cell_id=11, fill=dense_fuel, region=-r0) - c11.temperature = [500, 0, 700, 800] - c12 = openmc.Cell(cell_id=12, fill=moderator, region=+r0) - fuel_univ = openmc.Universe(universe_id=11, cells=(c11, c12)) - - lat = openmc.RectLattice(lattice_id=101) - lat.dimension = [2, 2] - lat.lower_left = [-2.0, -2.0] - lat.pitch = [2.0, 2.0] - lat.universes = [[fuel_univ]*2]*2 - lat.outer = mod_univ - - x0 = openmc.XPlane(x0=-3.0) - x1 = openmc.XPlane(x0=3.0) - y0 = openmc.YPlane(y0=-3.0) - y1 = openmc.YPlane(y0=3.0) - for s in [x0, x1, y0, y1]: - s.boundary_type = 'reflective' - c101 = openmc.Cell(cell_id=101, fill=lat, region=+x0 & -x1 & +y0 & -y1) - root_univ = openmc.Universe(universe_id=0, cells=(c101,)) - - geometry = openmc.Geometry(root_univ) - geometry.export_to_xml() - - #################### - # Settings - #################### - - sets_file = openmc.Settings() - sets_file.batches = 5 - sets_file.inactive = 0 - sets_file.particles = 1000 - sets_file.source = openmc.Source(space=openmc.stats.Box( - [-1, -1, -1], [1, 1, 1])) - sets_file.temperature = {'tolerance': 1000, 'multipole': True} - sets_file.export_to_xml() - - #################### - # Plots - #################### - - plot1 = openmc.Plot(plot_id=1) - plot1.basis = 'xy' - plot1.color_by = 'cell' - plot1.filename = 'cellplot' - plot1.origin = (0, 0, 0) - plot1.width = (7, 7) - plot1.pixels = (400, 400) - - plot2 = openmc.Plot(plot_id=2) - plot2.basis = 'xy' - plot2.color_by = 'material' - plot2.filename = 'matplot' - plot2.origin = (0, 0, 0) - plot2.width = (7, 7) - plot2.pixels = (400, 400) - - plots = openmc.Plots([plot1, plot2]) - plots.export_to_xml() - def execute_test(self): if not 'OPENMC_MULTIPOLE_LIBRARY' in os.environ: raise RuntimeError("The 'OPENMC_MULTIPOLE_LIBRARY' environment " @@ -107,5 +82,6 @@ class MultipoleTestHarness(PyAPITestHarness): if __name__ == '__main__': - harness = MultipoleTestHarness('statepoint.5.h5') + model = make_model() + harness = MultipoleTestHarness('statepoint.5.h5', model) harness.main() From 7d3b049fba641d1a4de29d8684b31a215f5c7445 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Apr 2017 19:39:08 -0500 Subject: [PATCH 2/3] Assume long long support by default (don't rely on fancy CMake CXX_STANDARD) --- CMakeLists.txt | 6 ------ src/pugixml/pugiconfig.hpp | 2 +- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2dc30a21d4..19d2d3886b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -235,12 +235,6 @@ endif() #=============================================================================== add_library(pugixml src/pugixml/pugixml_c.cpp src/pugixml/pugixml.cpp) -if(CMAKE_VERSION VERSION_LESS 3.1) - target_compile_options(pugixml PRIVATE -std=c++11) -else() - set_property(TARGET pugixml PROPERTY CXX_STANDARD 11) -endif() - add_library(pugixml_fortran src/pugixml/pugixml_f.F90) target_link_libraries(pugixml_fortran pugixml) diff --git a/src/pugixml/pugiconfig.hpp b/src/pugixml/pugiconfig.hpp index 2d1e6aa8b0..2f1027aa38 100644 --- a/src/pugixml/pugiconfig.hpp +++ b/src/pugixml/pugiconfig.hpp @@ -44,7 +44,7 @@ // #define PUGIXML_HEADER_ONLY // Uncomment this to enable long long support -// #define PUGIXML_HAS_LONG_LONG +#define PUGIXML_HAS_LONG_LONG #endif From c6202eb781bec10d6af3685b3ab3ae3cc1f8603e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Apr 2017 09:39:05 -0500 Subject: [PATCH 3/3] Update inputs for test_diff_tally --- tests/test_diff_tally/inputs_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat index 5667969957..aaf96c768f 100644 --- a/tests/test_diff_tally/inputs_true.dat +++ b/tests/test_diff_tally/inputs_true.dat @@ -306,7 +306,7 @@ -160 -160 -183 160 160 183 - True + true