From 535d823ddf9b51a425a76f2567d1c731cc05496b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 29 Jun 2021 15:04:41 +0700 Subject: [PATCH 1/3] Add temperature interpolation test --- tests/unit_tests/test_temp_interp.py | 155 +++++++++++++++++++++++++++ 1 file changed, 155 insertions(+) create mode 100644 tests/unit_tests/test_temp_interp.py diff --git a/tests/unit_tests/test_temp_interp.py b/tests/unit_tests/test_temp_interp.py new file mode 100644 index 0000000000..9d5bbdebd1 --- /dev/null +++ b/tests/unit_tests/test_temp_interp.py @@ -0,0 +1,155 @@ +from math import isnan +import os +from pathlib import Path + +from uncertainties import ufloat +import numpy as np +import openmc.data +from openmc.data import K_BOLTZMANN +from openmc.stats import Uniform +import pytest + + +def make_fake_cross_section(): + """Create fake U235 nuclide + + This nuclide is designed to have k_inf=1 at 300 K, k_inf=2 at 600 K, and + k_inf=1 at 900 K. The absorption cross section is also constant with + temperature so as to make the true k-effective go linear with temperature. + """ + + def isotropic_angle(): + return openmc.data.AngleDistribution( + [E_min, E_max], + [Uniform(-1., 1.), Uniform(-1., 1.)] + ) + + def cross_section(value): + return openmc.data.Tabulated1D( + energy, + value*np.ones_like(energy) + ) + + temperatures = (300, 600, 900) + + u235_fake = openmc.data.IncidentNeutron( + 'U235', 92, 235, 0, 233.0248, [T*K_BOLTZMANN for T in temperatures] + ) + + # Create energy grids + E_min, E_max = 1e-5, 20.0e6 + energy = np.logspace(np.log10(E_min), np.log10(E_max)) + for T in temperatures: + u235_fake.energy[f'{T}K'] = energy + + # Create elastic scattering + elastic = openmc.data.Reaction(2) + for T in temperatures: + elastic.xs[f'{T}K'] = cross_section(1.0) + elastic_dist = openmc.data.UncorrelatedAngleEnergy(isotropic_angle()) + product = openmc.data.Product() + product.distribution.append(elastic_dist) + elastic.products.append(product) + u235_fake.reactions[2] = elastic + + # Create fission + fission = openmc.data.Reaction(18) + fission.center_of_mass = False + fission.Q_value = 193.0e6 + fission_xs = (2., 4., 2.) + for T, xs in zip(temperatures, fission_xs): + fission.xs[f'{T}K'] = cross_section(xs) + a = openmc.data.Tabulated1D([E_min, E_max], [0.988e6, 0.988e6]) + b = openmc.data.Tabulated1D([E_min, E_max], [2.249e-6, 2.249e-6]) + fission_dist = openmc.data.UncorrelatedAngleEnergy( + isotropic_angle(), + openmc.data.WattEnergy(a, b, -E_max) + ) + product = openmc.data.Product() + product.distribution.append(fission_dist) + product.yield_ = openmc.data.Polynomial((2.0,)) + fission.products.append(product) + u235_fake.reactions[18] = fission + + # Create capture + capture = openmc.data.Reaction(102) + capture.q_value = 6.5e6 + capture_xs = (2., 0., 2.) + for T, xs in zip(temperatures, capture_xs): + capture.xs[f'{T}K'] = cross_section(xs) + u235_fake.reactions[102] = capture + + # Export HDF5 file + u235_fake.export_to_hdf5('U235_fake.h5', 'w') + + lib = openmc.data.DataLibrary() + lib.register_file('U235_fake.h5') + lib.export_to_xml('cross_sections_fake.xml') + + +@pytest.fixture(scope='module') +def model(tmp_path_factory): + tmp_path = tmp_path_factory.mktemp("temp_interp") + orig = Path.cwd() + os.chdir(tmp_path) + + make_fake_cross_section() + + model = openmc.model.Model() + mat = openmc.Material() + mat.add_nuclide('U235', 1.0) + model.materials.append(mat) + model.materials.cross_sections = str(Path('cross_sections_fake.xml').resolve()) + + 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.inactive = 0 + model.settings.batches = 10 + + tally = openmc.Tally() + tally.scores = ['absorption', 'fission', 'scatter', 'nu-fission'] + model.tallies = [tally] + + try: + yield model + finally: + os.chdir(orig) + + +@pytest.mark.parametrize( + ["method", "temperature", "fission_expected"], + [ + ("nearest", 300.0, 0.5), + ("nearest", 600.0, 1.0), + ("nearest", 900.0, 0.5), + ("interpolation", 360.0, 0.6), + ("interpolation", 450.0, 0.75), + ("interpolation", 540.0, 0.9), + ("interpolation", 660.0, 0.9), + ("interpolation", 750.0, 0.75), + ("interpolation", 840.0, 0.6), + ] +) +def test_interpolation(model, method, temperature, fission_expected): + model.settings.temperature = {'method': method, 'default': temperature} + sp_filename = model.run() + with openmc.StatePoint(sp_filename) as sp: + t = sp.tallies[model.tallies[0].id] + absorption_mean, fission_mean, scatter_mean, nu_fission_mean = t.mean.ravel() + absorption_unc, fission_unc, scatter_unc, nu_fission_unc = t.std_dev.ravel() + + nu = 2.0 + assert abs(absorption_mean - 1) < 3*absorption_unc + assert abs(fission_mean - fission_expected) < 3*fission_unc + assert abs(scatter_mean - 1/4) < 3*scatter_unc + assert abs(nu_fission_mean - nu*fission_expected) < 3*nu_fission_unc + + # Check that k-effective value matches expected + k = sp.k_combined + if isnan(k.s): + assert k.n == pytest.approx(nu*fission_expected) + else: + assert abs(k.n - nu*fission_expected) <= 3*k.s From 763c3a151ea7bb7cc7a8e08876db456006e82071 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Jul 2021 12:57:30 -0500 Subject: [PATCH 2/3] Update neighbor list documentation. Closes #1833 --- docs/source/methods/geometry.rst | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 04e3456c4d..2c7b7ecfb2 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -637,16 +637,20 @@ unsuccessful, then a search is done over every cell in the base universe. Building Neighbor Lists ----------------------- -After the geometry has been loaded and stored in memory from an input file, -OpenMC builds a list for each surface containing any cells that are bounded by -that surface in order to speed up processing of surface crossings. The algorithm -to build these lists is as follows. First, we loop over all cells in the -geometry and count up how many times each surface appears in a specification as -bounding a negative half-space and bounding a positive half-space. Two arrays -are then allocated for each surface, one that lists each cell that contains the -negative half-space of the surface and one that lists each cell that contains -the positive half-space of the surface. Another loop is performed over all cells -and the neighbor lists are populated for each surface. +Neighbor lists are data structures that are used to accelerate geometry searches +when a particle crosses a boundary. Namely, they are used to constrain the +number of cells that must be searched in order to determine which cell a +particle is crossing into. Earlier versions of OpenMC relied on "surface-based" +neighbor lists, where the cells that are adjacent to each surface are stored in +lists, one for each side of a surface. As of version 0.11, OpenMC switched to +using "cell-based" neighbor lists. For each cell, a list of the adjacent cells +is stored and then used to limit future searches. Unlike surface-based neighbor +lists, cell-based neighbor lists cannot be computed prior to transport. Thus, +cell-based neighbor lists in OpenMC grow dynamically as particles are +transported through the geometry and cross surfaces. Special care must be taken +to ensure that these dynamic neighbor lists are populated in a threadsafe +manner. Full details of the implementation in OpenMC can be found in a paper by +`Harper et al `_. .. _reflection: From 48be9c1c338d78ebf8b3b7345d4064ee32146a10 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Jul 2021 10:53:39 -0500 Subject: [PATCH 3/3] Respond to @pshriwise comments on #1858 --- tests/unit_tests/test_temp_interp.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/tests/unit_tests/test_temp_interp.py b/tests/unit_tests/test_temp_interp.py index 9d5bbdebd1..6cbe2f63fb 100644 --- a/tests/unit_tests/test_temp_interp.py +++ b/tests/unit_tests/test_temp_interp.py @@ -2,7 +2,6 @@ from math import isnan import os from pathlib import Path -from uncertainties import ufloat import numpy as np import openmc.data from openmc.data import K_BOLTZMANN @@ -18,7 +17,7 @@ def make_fake_cross_section(): temperature so as to make the true k-effective go linear with temperature. """ - def isotropic_angle(): + def isotropic_angle(E_min, E_max): return openmc.data.AngleDistribution( [E_min, E_max], [Uniform(-1., 1.), Uniform(-1., 1.)] @@ -40,13 +39,13 @@ def make_fake_cross_section(): E_min, E_max = 1e-5, 20.0e6 energy = np.logspace(np.log10(E_min), np.log10(E_max)) for T in temperatures: - u235_fake.energy[f'{T}K'] = energy + u235_fake.energy['{}K'.format(T)] = energy # Create elastic scattering elastic = openmc.data.Reaction(2) for T in temperatures: - elastic.xs[f'{T}K'] = cross_section(1.0) - elastic_dist = openmc.data.UncorrelatedAngleEnergy(isotropic_angle()) + elastic.xs['{}K'.format(T)] = cross_section(1.0) + elastic_dist = openmc.data.UncorrelatedAngleEnergy(isotropic_angle(E_min, E_max)) product = openmc.data.Product() product.distribution.append(elastic_dist) elastic.products.append(product) @@ -58,11 +57,11 @@ def make_fake_cross_section(): fission.Q_value = 193.0e6 fission_xs = (2., 4., 2.) for T, xs in zip(temperatures, fission_xs): - fission.xs[f'{T}K'] = cross_section(xs) + fission.xs['{}K'.format(T)] = cross_section(xs) a = openmc.data.Tabulated1D([E_min, E_max], [0.988e6, 0.988e6]) b = openmc.data.Tabulated1D([E_min, E_max], [2.249e-6, 2.249e-6]) fission_dist = openmc.data.UncorrelatedAngleEnergy( - isotropic_angle(), + isotropic_angle(E_min, E_max), openmc.data.WattEnergy(a, b, -E_max) ) product = openmc.data.Product() @@ -76,7 +75,7 @@ def make_fake_cross_section(): capture.q_value = 6.5e6 capture_xs = (2., 0., 2.) for T, xs in zip(temperatures, capture_xs): - capture.xs[f'{T}K'] = cross_section(xs) + capture.xs['{}K'.format(T)] = cross_section(xs) u235_fake.reactions[102] = capture # Export HDF5 file