crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'openmc-dev:develop' into parametric_tokamak_source

Browse source
This commit is contained in:
Ethan Peterson 2026-07-06 13:35:54 -04:00 • committed by GitHub
commit 23572ccc92
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
18 changed files with 553 additions and 64 deletions
Download patch file Download diff file Expand all files Collapse all files

12
include/openmc/math_functions.h
View file

@ -247,5 +247,17 @@ void get_energy_index(
double standard_normal_cdf(double z);
//==============================================================================
//! Return true if two floating-point values are approximately equal within a
//! combined relative and absolute tolerance.
//!
//! \param a first floating point value
//! \param b second floating point value
//! \param rel_tol relative tolerance
//! \param abs_tol absolute tolerance
//! \return true if a and b are approximately equal, false otherwise
//==============================================================================
bool isclose(double a, double b, double rel_tol, double abs_tol);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H

8
include/openmc/weight_windows.h
View file

@ -52,8 +52,12 @@ struct WeightWindow {
double weight_cutoff {DEFAULT_WEIGHT_CUTOFF};
int max_split {10};
//! Whether the weight window is in a valid state
bool is_valid() const { return lower_weight >= 0.0; }
//! Whether the weight window is in a valid state. A non-positive lower
//! bound indicates that no weight window information exists at this
//! location (generators mark such cells with -1, and a lower bound of zero
//! conventionally turns the weight window game off in a cell, as in MCNP
//! wwinp files), in which case no weight window game is played.
bool is_valid() const { return lower_weight > 0.0; }
//! Adjust the weight window by a constant factor
void scale(double factor)

129
openmc/data/dose/dose.py
View file

@ -4,50 +4,76 @@ import numpy as np
import openmc.checkvalue as cv
_FILES = {
('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt',
('icrp74', 'photon'): Path('icrp74') / 'photons.txt',
('icrp116', 'electron'): Path('icrp116') / 'electrons.txt',
('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt',
('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt',
('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt',
('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt',
('icrp116', 'photon'): Path('icrp116') / 'photons.txt',
('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt',
('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt',
('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt',
('icrp116', 'positron'): Path('icrp116') / 'positrons.txt',
('icrp116', 'proton'): Path('icrp116') / 'protons.txt',
_FULL_GEOMETRIES = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
_LIMITED_GEOMETRIES = ('AP', 'PA', 'ISO')
_TABLES = {
('icrp74', 'effective', 'neutron'): (
Path('icrp74') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp74', 'effective', 'photon'): (
Path('icrp74') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'electron'): (
Path('icrp116') / 'electrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'helium'): (
Path('icrp116') / 'helium_ions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'mu-'): (
Path('icrp116') / 'negative_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi-'): (
Path('icrp116') / 'negative_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'neutron'): (
Path('icrp116') / 'neutrons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon'): (
Path('icrp116') / 'photons.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'photon kerma'): (
Path('icrp116') / 'photons_kerma.txt', _FULL_GEOMETRIES),
('icrp116', 'effective', 'mu+'): (
Path('icrp116') / 'positive_muons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'pi+'): (
Path('icrp116') / 'positive_pions.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'positron'): (
Path('icrp116') / 'positrons.txt', _LIMITED_GEOMETRIES),
('icrp116', 'effective', 'proton'): (
Path('icrp116') / 'protons.txt', _FULL_GEOMETRIES),
('icrp74', 'ambient', 'neutron'): (
Path('icrp74') / 'neutrons_H10.txt', None),
('icrp74', 'ambient', 'photon'): (
Path('icrp74') / 'photons_H10.txt', None),
}
_DOSE_TABLES = {}
def _load_dose_icrp(data_source: str, particle: str):
"""Load effective dose tables from text files.
def _load_dose_table(data_source: str, dose_quantity: str, particle: str):
"""Load dose tables from text files.
Parameters
----------
data_source : {'icrp74', 'icrp116'}
The dose conversion data source to use
dose_quantity : {'effective', 'ambient'}
Dose quantity to load. 'ambient' corresponds to ambient dose
equivalent H*(10).
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
"""
path = Path(__file__).parent / _FILES[data_source, particle]
key = (data_source, dose_quantity, particle)
path = Path(__file__).parent / _TABLES[key][0]
data = np.loadtxt(path, skiprows=3, encoding='utf-8')
data[:, 0] *= 1e6 # Change energies to eV
_DOSE_TABLES[data_source, particle] = data
_DOSE_TABLES[key] = data
def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
"""Return effective dose conversion coefficients.
def dose_coefficients(
particle, geometry='AP', data_source='icrp116', dose_quantity='effective'
):
"""Return dose conversion coefficients.
This function provides fluence (and air kerma) to effective or ambient dose
(H*(10)) conversion coefficients for various types of external exposures
based on values in ICRP publications. Corrected values found in a
corrigendum are used rather than the values in the original report.
Available libraries include `ICRP Publication 74
This function provides fluence (and air kerma) to effective dose or ambient
dose equivalent (H*(10)) conversion coefficients for various types of
external exposures based on values in ICRP publications. Corrected values
found in a corrigendum are used rather than the values in the original
report. Available libraries include `ICRP Publication 74
<https://doi.org/10.1016/S0146-6453(96)90010-X>` and `ICRP Publication 116
<https://doi.org/10.1016/j.icrp.2011.10.001>`.
@ -63,45 +89,58 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'):
particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'}
Incident particle
geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}
Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the
meaning of the options here.
Irradiation geometry assumed for effective dose coefficients. Refer to
ICRP-116 (Section 3.2) for the meaning of the options here. This
argument does not apply when ``dose_quantity`` is 'ambient'.
data_source : {'icrp74', 'icrp116'}
The data source for the effective dose conversion coefficients.
The data source for the dose conversion coefficients.
dose_quantity : {'effective', 'ambient'}
Dose quantity to return. 'effective' returns effective dose
coefficients; 'ambient' returns ambient dose equivalent (H*(10))
coefficients.
Returns
-------
energy : numpy.ndarray
Energies at which dose conversion coefficients are given
dose_coeffs : numpy.ndarray
Effective dose coefficients in [pSv cm^2] at provided energies. For
'photon kerma', the coefficients are given in [Sv/Gy].
Dose coefficients in [pSv cm^2] at provided energies. For 'photon
kerma', the coefficients are given in [Sv/Gy].
"""
cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'})
cv.check_value('geometry', geometry, _FULL_GEOMETRIES)
cv.check_value('data_source', data_source, {'icrp74', 'icrp116'})
cv.check_value('dose_quantity', dose_quantity, {'effective', 'ambient'})
if (data_source, particle) not in _FILES:
available_particles = sorted({p for (ds, p) in _FILES if ds == data_source})
key = (data_source, dose_quantity, particle)
if key not in _TABLES:
available_particles = sorted(
p for ds, dq, p in _TABLES
if ds == data_source and dq == dose_quantity
)
msg = (
f"'{particle}' has no dose data in data source {data_source}. "
f"Available particles for {data_source} are: {available_particles}"
f"'{particle}' has no {dose_quantity} dose data in data source "
f"{data_source}. Available particles for {data_source} "
f"with dose quantity {dose_quantity} are: {available_particles}"
)
raise ValueError(msg)
elif (data_source, particle) not in _DOSE_TABLES:
_load_dose_icrp(data_source, particle)
elif key not in _DOSE_TABLES:
_load_dose_table(data_source, dose_quantity, particle)
# Get all data for selected particle
data = _DOSE_TABLES[data_source, particle]
data = _DOSE_TABLES[key]
columns = _TABLES[key][1]
# Determine index for selected geometry
if particle in ('neutron', 'photon', 'proton', 'photon kerma'):
columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO')
if columns is None:
if geometry != 'AP':
raise ValueError(
"Irradiation geometry is not defined for ambient dose "
"equivalent coefficients. Use the default geometry='AP'."
)
index = 0
else:
columns = ('AP', 'PA', 'ISO')
index = columns.index(geometry)
index = columns.index(geometry)
# Pull out energy and dose from table
energy = data[:, 0].copy()
dose_coeffs = data[:, index + 1].copy()
return energy, dose_coeffs

50
openmc/data/dose/icrp74/neutrons_H10.txt Normal file
View file

@ -0,0 +1,50 @@
Neutrons: Ambient per fluence, in units of pSv cm², for monoenergetic particles incident.
Energy (MeV) Dose
1.00E-09 6.60
1.00E-08 9.00
2.53E-08 10.6
1.00E-07 12.9
2.00E-07 13.5
5.00E-07 13.6
1.00E-06 13.3
2.00E-06 12.9
5.00E-06 12.0
1.00E-05 11.3
2.00E-05 10.6
5.00E-05 9.90
1.00E-04 9.40
2.00E-04 8.90
5.00E-04 8.30
1.00E-03 7.90
2.00E-03 7.70
5.00E-03 8.00
1.00E-02 10.5
2.00E-02 16.6
3.00E-02 23.7
5.00E-02 41.1
7.00E-02 60.0
1.00E-01 88.0
1.50E-01 132
2.00E-01 170
3.00E-01 233
5.00E-01 322
7.00E-01 375
9.00E-01 400
1 416
1.2 425
2 420
3 412
4 408
5 405
6 400
7 405
8 409
9 420
10 440
12 480
14 520
15 540
16 555
18 570
20 600

28
openmc/data/dose/icrp74/photons_H10.txt Normal file
View file

@ -0,0 +1,28 @@
Photons: Ambient dose (H*10) per fluence, in units of pSv cm²
Energy (MeV) Dose
0.010 0.061
0.015 0.83
0.020 1.05
0.030 0.81
0.040 0.64
0.050 0.55
0.060 0.51
0.080 0.53
0.100 0.61
0.150 0.89
0.200 1.20
0.300 1.80
0.400 2.38
0.500 2.93
0.600 3.44
0.800 4.38
1 5.20
1.5 6.90
2 8.60
3 11.1
4 13.4
5 15.5
6 17.6
8 21.6
10 25.6

3
openmc/deplete/microxs.py
View file

@ -233,7 +233,8 @@ def get_microxs_and_flux(
model.export_to_model_xml()
comm.barrier()
# Reinitialize with tallies
openmc.lib.init(intracomm=comm)
output = run_kwargs.get('output', True) if run_kwargs else True
openmc.lib.init(intracomm=comm, output=output)
with TemporaryDirectory() as temp_dir:
# Indicate to run in temporary directory unless being executed through

2
openmc/deplete/r2s.py
View file

@ -347,7 +347,7 @@ class R2SManager:
# Run neutron transport and get fluxes and micros. Run via openmc.lib to
# maintain a consistent parallelism strategy with the activation step.
with TemporarySession():
with TemporarySession(output=False):
self.results['fluxes'], self.results['micros'] = get_microxs_and_flux(
self.neutron_model, domains, **micro_kwargs)

33
src/lattice.cpp
View file

@ -10,6 +10,7 @@
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/string_utils.h"
#include "openmc/vector.h"
#include "openmc/xml_interface.h"
@ -260,25 +261,33 @@ std::pair<double, array<int, 3>> RectLattice::distance(
// Determine the oncoming edge.
double x0 {copysign(0.5 * pitch_[0], u.x)};
double y0 {copysign(0.5 * pitch_[1], u.y)};
double z0;
double d = std::min(
u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY);
// Evaluate the distance to each oncoming edge independently. Comparing these
// distances directly (rather than reconstructing the crossing position)
// avoids the floating-point cancellation that occurs for large pitches.
double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY;
double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY;
double dz = INFTY;
if (is_3d_) {
z0 = copysign(0.5 * pitch_[2], u.z);
d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY);
double z0 {copysign(0.5 * pitch_[2], u.z)};
dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY;
}
// Determine which lattice boundaries are being crossed
// The distance to the nearest lattice boundary is the smallest axial
// distance.
double d = std::min({dx, dy, dz});
// Determine which lattice boundaries are being crossed. The axis attaining
// the minimum is exactly equal to d, so at least one translation is always
// set for a finite crossing; a near-equal second axis indicates a corner
// crossing.
array<int, 3> lattice_trans = {0, 0, 0};
if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION)
if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION))
lattice_trans[0] = copysign(1, u.x);
if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION)
if (isclose(d, dy, FP_COINCIDENT, FP_PRECISION))
lattice_trans[1] = copysign(1, u.y);
if (is_3d_) {
if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION)
lattice_trans[2] = copysign(1, u.z);
}
if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION))
lattice_trans[2] = copysign(1, u.z);
return {d, lattice_trans};
}

8
src/math_functions.cpp
View file

@ -1001,4 +1001,12 @@ void get_energy_index(
}
}
// Return true if two floating-point values are approximately equal within a
// combined relative and absolute tolerance.
bool isclose(double a, double b, double rel_tol, double abs_tol)
{
return std::abs(a - b) <=
std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol);
}
} // namespace openmc

3
src/particle.cpp
View file

@ -338,13 +338,14 @@ void Particle::event_cross_surface()
boundary().lattice_translation()[2] != 0) {
// Particle crosses lattice boundary
int i_lattice = coord(boundary().coord_level() - 1).lattice();
bool verbose = settings::verbosity >= 10 || trace();
cross_lattice(*this, boundary(), verbose);
event() = TallyEvent::LATTICE;
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
auto& lat {*model::lattices[lowest_coord().lattice()]};
auto& lat {*model::lattices[i_lattice]};
bool is_valid;
Direction normal =
lat.get_normal(boundary().lattice_translation(), is_valid);

21
tests/cpp_unit_tests/test_math.cpp
View file

@ -355,3 +355,24 @@ TEST_CASE("Test broaden_wmp_polynomials")
REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val));
}
}
TEST_CASE("Test isclose")
{
using openmc::isclose;
// Identical values are always close, regardless of tolerances.
REQUIRE(isclose(1.0, 1.0, 0.0, 0.0));
REQUIRE(isclose(0.0, 0.0, 0.0, 0.0));
// Absolute tolerance governs comparisons near zero.
REQUIRE(isclose(0.0, 1e-15, 0.0, 1e-14));
REQUIRE_FALSE(isclose(0.0, 1e-13, 0.0, 1e-14));
// Relative tolerance scales with the magnitude of the operands.
REQUIRE(isclose(1.0e12, 1.0e12 + 1.0, 1e-12, 0.0));
REQUIRE_FALSE(isclose(1.0e12, 1.0e12 + 10.0, 1e-12, 0.0));
// The looser of the two tolerances wins.
REQUIRE(isclose(1.0, 1.0 + 1e-13, 0.0, 1e-12));
REQUIRE(isclose(1.0e6, 1.0e6 + 1e-4, 1e-9, 0.0));
}

0
tests/regression_tests/lattice_large_pitch/__init__.py Normal file
View file

51
tests/regression_tests/lattice_large_pitch/inputs_true.dat Normal file
View file

@ -0,0 +1,51 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<material id="1">
<density value="0.001" units="g/cm3"/>
<nuclide name="N14" ao="1.0"/>
</material>
<material id="2">
<density value="7.0" units="g/cm3"/>
<nuclide name="Fe56" ao="1.0"/>
</material>
</materials>
<geometry>
<cell id="1" material="2" universe="1"/>
<cell id="2" material="1" universe="2"/>
<cell id="3" fill="3" region="1 -2 3 -4" universe="4"/>
<lattice id="3">
<pitch>100.0 100.0</pitch>
<outer>2</outer>
<dimension>3 3</dimension>
<lower_left>-150.0 -150.0</lower_left>
<universes>
1 2 1
2 1 2
1 2 1 </universes>
</lattice>
<surface id="1" name="minimum x" type="x-plane" boundary="vacuum" coeffs="-150.0"/>
<surface id="2" name="maximum x" type="x-plane" boundary="vacuum" coeffs="150.0"/>
<surface id="3" name="minimum y" type="y-plane" boundary="vacuum" coeffs="-150.0"/>
<surface id="4" name="maximum y" type="y-plane" boundary="vacuum" coeffs="150.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>1000</particles>
<batches>10</batches>
</settings>
<tallies>
<mesh id="1">
<dimension>6 6</dimension>
<lower_left>-150.0 -150.0</lower_left>
<upper_right>150.0 150.0</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

73
tests/regression_tests/lattice_large_pitch/results_true.dat Normal file
View file

@ -0,0 +1,73 @@
tally 1:
1.749265E+01
3.275925E+01
4.159430E+01
1.799568E+02
7.027454E+01
4.967568E+02
6.789716E+01
4.652258E+02
4.272935E+01
1.861129E+02
2.073342E+01
4.374279E+01
4.028297E+01
1.641642E+02
6.734263E+01
4.627171E+02
1.037657E+02
1.078682E+03
1.108394E+02
1.240927E+03
7.236770E+01
5.302034E+02
4.356361E+01
1.921990E+02
6.731536E+01
4.615977E+02
9.850684E+01
9.810272E+02
5.164863E+02
2.670336E+04
5.097770E+02
2.604770E+04
1.064847E+02
1.137536E+03
7.052986E+01
5.003220E+02
7.065046E+01
5.036723E+02
1.044890E+02
1.103250E+03
5.121544E+02
2.632389E+04
5.065331E+02
2.570752E+04
1.044794E+02
1.101249E+03
6.569142E+01
4.361120E+02
4.739812E+01
2.313289E+02
7.402284E+01
5.591689E+02
1.066905E+02
1.149521E+03
1.038665E+02
1.086813E+03
7.160008E+01
5.217744E+02
4.219705E+01
1.816100E+02
2.089838E+01
4.464899E+01
4.154271E+01
1.745866E+02
7.081253E+01
5.049997E+02
6.673273E+01
4.509038E+02
4.184624E+01
1.771546E+02
1.853898E+01
3.538608E+01

70
tests/regression_tests/lattice_large_pitch/test.py Normal file
View file

@ -0,0 +1,70 @@
"""Regression test for rectangular lattices with large pitch values.
Large pitches used to trigger a segmentation fault in ``RectLattice::distance``
because the boundary-crossing check compared a reconstructed crossing position
against an absolute tolerance that did not scale with the geometry (see #3852).
This test transports particles across a lattice with a large pitch to ensure the
crossing logic remains robust.
"""
import openmc
import pytest
from tests.testing_harness import PyAPITestHarness
@pytest.fixture
def model():
model = openmc.Model()
# Large pitch that previously caused floating-point cancellation
pitch = 100.0
n = 3
air = openmc.Material()
air.set_density('g/cm3', 0.001)
air.add_nuclide('N14', 1.0)
metal = openmc.Material()
metal.set_density('g/cm3', 7.0)
metal.add_nuclide('Fe56', 1.0)
metal_cell = openmc.Cell(fill=metal)
metal_uni = openmc.Universe(cells=[metal_cell])
air_cell = openmc.Cell(fill=air)
air_uni = openmc.Universe(cells=[air_cell])
lattice = openmc.RectLattice()
lattice.lower_left = (-pitch*n/2, -pitch*n/2)
lattice.pitch = (pitch, pitch)
lattice.outer = air_uni
lattice.universes = [
[metal_uni, air_uni, metal_uni],
[air_uni, metal_uni, air_uni],
[metal_uni, air_uni, metal_uni],
]
box = openmc.model.RectangularPrism(pitch*n, pitch*n, boundary_type='vacuum')
root_cell = openmc.Cell(region=-box, fill=lattice)
model.geometry = openmc.Geometry([root_cell])
model.settings.run_mode = 'fixed source'
model.settings.batches = 10
model.settings.particles = 1000
mesh = openmc.RegularMesh()
mesh.dimension = (6, 6)
mesh.lower_left = (-pitch*n/2, -pitch*n/2)
mesh.upper_right = (pitch*n/2, pitch*n/2)
mesh_filter = openmc.MeshFilter(mesh)
tally = openmc.Tally(tally_id=1)
tally.filters = [mesh_filter]
tally.scores = ['flux']
model.tallies = [tally]
return model
def test_lattice_large_pitch(model):
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

55
tests/regression_tests/weightwindows/test.py
View file

@ -119,6 +119,61 @@ def test_weightwindows(shared_secondary, subdir):
test.main()
def test_zero_bound_windows_play_no_game(tmp_path):
# A weight window lower bound of zero means no weight window information
# exists there (MCNP wwinp files use zero to turn the game off in a cell),
# so transport must proceed as if weight windows were disabled. Previously,
# zero-bound windows demanded a split at every checkpoint (weight/0 ->
# max_split), multiplying the particle population until terminated by the
# split or weight cutoff limits.
model = build_model(False)
for ww in model.settings.weight_windows:
ww.lower_ww_bounds = np.zeros_like(ww.lower_ww_bounds)
ww.upper_ww_bounds = np.zeros_like(ww.upper_ww_bounds)
sp_zero = model.run(cwd=tmp_path / 'zero_windows')
model.settings.weight_windows_on = False
sp_off = model.run(cwd=tmp_path / 'windows_off')
with openmc.StatePoint(sp_zero) as sp:
flux_zero = list(sp.tallies.values())[0].mean
with openmc.StatePoint(sp_off) as sp:
flux_off = list(sp.tallies.values())[0].mean
np.testing.assert_allclose(flux_zero, flux_off, rtol=1e-12)
def test_zero_and_negative_bounds_equivalent(tmp_path):
# Zero and negative lower bounds both mean that no weight window
# information exists in a cell (generators mark such cells with -1, and
# MCNP wwinp files use zero), so they must produce identical transport.
# Unlike the all-zero case above, here particles are born under valid
# windows and encounter the no-information region in flight; previously a
# zero lower bound in that situation demanded a split at every checkpoint
# in the cell (weight/0 -> max_split), multiplying the particle population,
# while -1 played no game.
def run_with(bound_value, subdir):
model = build_model(False)
for ww in model.settings.weight_windows:
lb = np.array(ww.lower_ww_bounds, copy=True)
ub = np.array(ww.upper_ww_bounds, copy=True)
lb[3:, :, :, :] = bound_value
ub[3:, :, :, :] = bound_value
ww.lower_ww_bounds = lb
ww.upper_ww_bounds = ub
return model.run(cwd=tmp_path / subdir)
sp_zero = run_with(0.0, 'zero_region')
sp_negative = run_with(-1.0, 'negative_region')
with openmc.StatePoint(sp_zero) as sp:
flux_zero = list(sp.tallies.values())[0].mean
with openmc.StatePoint(sp_negative) as sp:
flux_negative = list(sp.tallies.values())[0].mean
np.testing.assert_allclose(flux_zero, flux_negative, rtol=1e-12)
def test_wwinp_cylindrical():
ww = openmc.WeightWindowsList.from_wwinp('ww_n_cyl.txt')[0]

27
tests/unit_tests/test_data_dose.py
View file

@ -34,6 +34,20 @@ def test_dose_coefficients():
assert energy[-1] == approx(20.0e6)
assert dose[-1] == approx(338.0)
energy, dose = dose_coefficients(
'neutron', data_source='icrp74', dose_quantity='ambient')
assert energy[0] == approx(1e-3)
assert dose[0] == approx(6.60)
assert energy[-1] == approx(20.0e6)
assert dose[-1] == approx(600)
energy, dose = dose_coefficients(
'photon', data_source='icrp74', dose_quantity='ambient')
assert energy[0] == approx(0.01e6)
assert dose[0] == approx(0.061)
assert energy[-1] == approx(10e6)
assert dose[-1] == approx(25.6)
# Invalid particle/geometry should raise an exception
with raises(ValueError):
dose_coefficients('slime', 'LAT')
@ -41,6 +55,14 @@ def test_dose_coefficients():
dose_coefficients('neutron', 'ZZ')
with raises(ValueError):
dose_coefficients('neutron', data_source='icrp7000')
with raises(ValueError):
dose_coefficients('neutron', dose_quantity='banana')
with raises(ValueError):
dose_coefficients(
'neutron', data_source='icrp116', dose_quantity='ambient')
with raises(ValueError):
dose_coefficients(
'neutron', 'ISO', data_source='icrp74', dose_quantity='ambient')
with raises(ValueError) as excinfo:
dose_coefficients("photons", data_source="icrp116")
expected_particles = [
@ -57,7 +79,8 @@ def test_dose_coefficients():
"proton",
]
expected_msg = (
"'photons' has no dose data in data source icrp116. "
f"Available particles for icrp116 are: {expected_particles}"
"'photons' has no effective dose data in data source icrp116. "
"Available particles for icrp116 with dose quantity effective are: "
f"{expected_particles}"
)
assert str(excinfo.value) == expected_msg

44
tests/unit_tests/test_surface_flux.py
View file

@ -98,6 +98,50 @@ def test_surface_filter_flux_angled(two_cell_model, run_in_tmpdir):
assert flux_mean == pytest.approx(1.0 / mu)
def test_surface_tally_during_lattice_crossing(run_in_tmpdir):
openmc.reset_auto_ids()
model = openmc.Model()
xmin = openmc.XPlane(-1.0, boundary_type="vacuum")
xmax = openmc.XPlane(1.0, boundary_type="vacuum")
ymin = openmc.YPlane(-1.0, boundary_type="vacuum")
ymax = openmc.YPlane(1.0, boundary_type="vacuum")
zmin = openmc.ZPlane(-1.0, boundary_type="vacuum")
zmax = openmc.ZPlane(1.0, boundary_type="vacuum")
inner_cell = openmc.Cell()
inner_univ = openmc.Universe(cells=[inner_cell])
tile_cell = openmc.Cell(fill=inner_univ)
tile_univ = openmc.Universe(cells=[tile_cell])
lattice = openmc.RectLattice()
lattice.lower_left = (-1.0, -1.0)
lattice.pitch = (1.0, 2.0)
lattice.universes = [[tile_univ, tile_univ]]
root_cell = openmc.Cell(
fill=lattice, region=+xmin & -xmax & +ymin & -ymax & +zmin & -zmax)
model.geometry = openmc.Geometry([root_cell])
src = openmc.IndependentSource()
src.space = openmc.stats.Point((-0.5, 0.0, 0.0))
src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
model.settings.run_mode = 'fixed source'
model.settings.batches = 1
model.settings.particles = 5
model.settings.source = src
current_tally = openmc.Tally()
current_tally.filters = [openmc.SurfaceFilter(xmax)]
current_tally.scores = ['current']
model.tallies = [current_tally]
model.run(apply_tally_results=True)
assert current_tally.mean.flat[0] == pytest.approx(1.0)
def test_cellfrom_filter_flux_directional(two_cell_model, run_in_tmpdir):
"""SurfaceFilter + CellFromFilter + flux scores only the correct direction."""
model, xmid, cell1, cell2 = two_cell_model