From 953b35e742d4526b4bf4c5429d876f21506ec5c3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 07:21:08 -0500 Subject: [PATCH 1/5] Add openmc_sample_external_source function and Python binding --- openmc/lib/core.py | 39 ++++++++++++++++++++++++++++++++++++++- openmc/lib/math.py | 11 +++++------ src/source.cpp | 10 ++++++++++ 3 files changed, 53 insertions(+), 7 deletions(-) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index de5f4adf21..3f4a301cd8 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -1,8 +1,10 @@ from contextlib import contextmanager from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p, - c_char, POINTER, Structure, c_void_p, create_string_buffer) + c_char, POINTER, Structure, c_void_p, create_string_buffer, + c_uint64) import sys import os +from random import getrandbits import numpy as np from numpy.ctypeslib import as_array @@ -10,6 +12,7 @@ from numpy.ctypeslib import as_array from . import _dll from .error import _error_handler import openmc.lib +import openmc class _SourceSite(Structure): @@ -95,6 +98,9 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), POINTER(c_double)] _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler +_dll.openmc_sample_external_source.argtypes = [POINTER(c_uint64), POINTER(_SourceSite)] +_dll.openmc_sample_external_source.restype = c_int +_dll.openmc_sample_external_source.errcheck = _error_handler def global_bounding_box(): @@ -415,6 +421,37 @@ def run(output=True): _dll.openmc_run() +def sample_external_source(prn_seed=None): + """Sample external source + + .. versionadded:: 0.13.1 + + Parameters + ---------- + prn_seed : int + PRNG seed; if None, one will be generated randomly + + Returns + ------- + openmc.SourceParticle + Sampled source particle + + """ + if prn_seed is None: + prn_seed = getrandbits(63) + + # Call into C API to sample source + site = _SourceSite() + _dll.openmc_sample_external_source(c_uint64(prn_seed), site) + + # Convert to SourceParticle and return + return openmc.SourceParticle( + r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, + delayed_group=site.delayed_group, surf_id=site.surf_id, + particle=openmc.ParticleType(site.particle) + ) + + def simulation_init(): """Initialize simulation""" _dll.openmc_simulation_init() diff --git a/openmc/lib/math.py b/openmc/lib/math.py index bdb09b2d1f..d70c20a9f6 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -1,12 +1,11 @@ from ctypes import c_int, c_double, POINTER, c_uint64 +from random import getrandbits import numpy as np from numpy.ctypeslib import ndpointer from . import _dll -from random import getrandbits - _dll.t_percentile.restype = c_double _dll.t_percentile.argtypes = [c_double, c_int] @@ -240,10 +239,10 @@ def maxwell_spectrum(T, prn_seed=None): Sampled outgoing energy """ - + if prn_seed is None: prn_seed = getrandbits(63) - + return _dll.maxwell_spectrum(T, c_uint64(prn_seed)) @@ -265,7 +264,7 @@ def watt_spectrum(a, b, prn_seed=None): Sampled outgoing energy """ - + if prn_seed is None: prn_seed = getrandbits(63) @@ -290,7 +289,7 @@ def normal_variate(mean_value, std_dev, prn_seed=None): Sampled outgoing normally distributed value """ - + if prn_seed is None: prn_seed = getrandbits(63) diff --git a/src/source.cpp b/src/source.cpp index 11ede5fc4e..8b634ec459 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -399,4 +399,14 @@ void free_memory_source() model::external_sources.clear(); } +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_sample_external_source(uint64_t* seed, SourceSite* site) +{ + *site = sample_external_source(seed); + return 0; +} + } // namespace openmc From b1267f05448d80bce3069779f03de68749c44297 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 13:23:07 -0500 Subject: [PATCH 2/5] Allow openmc_sample_external_source to produce multiple samples --- docs/source/pythonapi/capi.rst | 1 + include/openmc/capi.h | 1 + openmc/lib/core.py | 33 ++++++++++++++++++++------------- src/source.cpp | 13 +++++++++++-- 4 files changed, 33 insertions(+), 15 deletions(-) diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 03be25abf4..882e3c71b3 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -29,6 +29,7 @@ Functions reset run run_in_memory + sample_external_source simulation_init simulation_finalize source_bank diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 0929a11f76..69d3ff1f61 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -121,6 +121,7 @@ int openmc_regular_mesh_set_params(int32_t index, int n, const double* ll, int openmc_reset(); int openmc_reset_timers(); int openmc_run(); +int openmc_sample_external_source(size_t n, uint64_t* seed, void* sites); void openmc_set_seed(int64_t new_seed); int openmc_set_n_batches( int32_t n_batches, bool set_max_batches, bool add_statepoint_batch); diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 3f4a301cd8..3a437a2952 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -1,7 +1,7 @@ from contextlib import contextmanager from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p, c_char, POINTER, Structure, c_void_p, create_string_buffer, - c_uint64) + c_uint64, c_size_t) import sys import os from random import getrandbits @@ -98,7 +98,7 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double), POINTER(c_double)] _dll.openmc_global_bounding_box.restype = c_int _dll.openmc_global_bounding_box.errcheck = _error_handler -_dll.openmc_sample_external_source.argtypes = [POINTER(c_uint64), POINTER(_SourceSite)] +_dll.openmc_sample_external_source.argtypes = [c_size_t, POINTER(c_uint64), POINTER(_SourceSite)] _dll.openmc_sample_external_source.restype = c_int _dll.openmc_sample_external_source.errcheck = _error_handler @@ -421,35 +421,42 @@ def run(output=True): _dll.openmc_run() -def sample_external_source(prn_seed=None): +def sample_external_source(n_samples=1, prn_seed=None): """Sample external source .. versionadded:: 0.13.1 Parameters ---------- + n_samples : int + Number of samples prn_seed : int PRNG seed; if None, one will be generated randomly Returns ------- - openmc.SourceParticle - Sampled source particle + list of openmc.SourceParticle + List of samples source particles """ + if n_samples <= 0: + raise ValueError("Number of samples must be positive") if prn_seed is None: prn_seed = getrandbits(63) # Call into C API to sample source - site = _SourceSite() - _dll.openmc_sample_external_source(c_uint64(prn_seed), site) + sites_array = (_SourceSite * n_samples)() + _dll.openmc_sample_external_source(c_size_t(n_samples), c_uint64(prn_seed), sites_array) - # Convert to SourceParticle and return - return openmc.SourceParticle( - r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, - delayed_group=site.delayed_group, surf_id=site.surf_id, - particle=openmc.ParticleType(site.particle) - ) + # Convert to list of SourceParticle and return + return [ + openmc.SourceParticle( + r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt, + delayed_group=site.delayed_group, surf_id=site.surf_id, + particle=openmc.ParticleType(site.particle) + ) + for site in sites_array + ] def simulation_init(): diff --git a/src/source.cpp b/src/source.cpp index 8b634ec459..3716720dff 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -403,9 +403,18 @@ void free_memory_source() // C API //============================================================================== -extern "C" int openmc_sample_external_source(uint64_t* seed, SourceSite* site) +extern "C" int openmc_sample_external_source( + size_t n, uint64_t* seed, void* sites) { - *site = sample_external_source(seed); + if (!sites || !seed) { + set_errmsg("Received null pointer."); + return OPENMC_E_INVALID_ARGUMENT; + } + + auto sites_array = static_cast(sites); + for (size_t i = 0; i < n; ++i) { + sites_array[i] = sample_external_source(seed); + } return 0; } From ba5f6f58fad0b28a2c69d4cc63d7dff0ac069067 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Jul 2022 19:57:58 -0500 Subject: [PATCH 3/5] Add test for openmc.lib.sample_external_source --- tests/unit_tests/test_lib.py | 41 ++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index a0f18f8513..f1f84bff22 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -807,3 +807,44 @@ def test_cell_rotation(pincell_model_w_univ, mpi_intracomm): cell.rotation = (180., 0., 0.) assert cell.rotation == pytest.approx([180., 0., 0.]) openmc.lib.finalize() + + +def test_sample_external_source(run_in_tmpdir, mpi_intracomm): + # Define a simple model and export + mat = openmc.Material() + mat.add_nuclide('U235', 1.0e-2) + sph = openmc.Sphere(r=100.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.Source( + space=openmc.stats.Box([-5., -5., -5.], [5., 5., 5.]), + angle=openmc.stats.Monodirectional((0., 0., 1.)), + energy=openmc.stats.Discrete([1.0e5], [1.0]) + ) + model.settings.particles = 1000 + model.settings.batches = 10 + model.export_to_xml() + + # Sample some particles and make sure they match specified source + openmc.lib.init() + particles = openmc.lib.sample_external_source(10, prn_seed=3) + for p in particles: + assert -5. < p.r[0] < 5. + assert -5. < p.r[1] < 5. + assert -5. < p.r[2] < 5. + assert p.u[0] == 0.0 + assert p.u[1] == 0.0 + assert p.u[2] == 1.0 + assert p.E == 1.0e5 + + # Using the same seed should produce the same particles + other_particles = openmc.lib.sample_external_source(10, prn_seed=3) + for p1, p2 in zip(particles, other_particles): + assert p1.r == p2.r + assert p1.u == p2.u + assert p1.E == p2.E + assert p1.time == p2.time + assert p1.wgt == p2.wgt + + openmc.lib.finalize() From d8425d2669b71efbda64a3fbdaf645cf8468e291 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jul 2022 06:48:51 -0500 Subject: [PATCH 4/5] Apply @shimwell suggestions from code review Co-authored-by: Jonathan Shimwell --- tests/unit_tests/test_lib.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index f1f84bff22..b572429835 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -829,6 +829,7 @@ def test_sample_external_source(run_in_tmpdir, mpi_intracomm): # Sample some particles and make sure they match specified source openmc.lib.init() particles = openmc.lib.sample_external_source(10, prn_seed=3) + assert len(particles) == 10 for p in particles: assert -5. < p.r[0] < 5. assert -5. < p.r[1] < 5. @@ -840,6 +841,7 @@ def test_sample_external_source(run_in_tmpdir, mpi_intracomm): # Using the same seed should produce the same particles other_particles = openmc.lib.sample_external_source(10, prn_seed=3) + assert len(other_particles) == 10 for p1, p2 in zip(particles, other_particles): assert p1.r == p2.r assert p1.u == p2.u From 476fdd1bf7169ac08d17cc529ae9fe1f8bbabcc8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 18 Jul 2022 06:49:14 -0500 Subject: [PATCH 5/5] Improve description of prn_seed argument in openmc.lib functions --- openmc/lib/core.py | 3 ++- openmc/lib/math.py | 12 ++++++++---- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 3a437a2952..479e414aeb 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -431,7 +431,8 @@ def sample_external_source(n_samples=1, prn_seed=None): n_samples : int Number of samples prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- diff --git a/openmc/lib/math.py b/openmc/lib/math.py index d70c20a9f6..0f61450bfd 100644 --- a/openmc/lib/math.py +++ b/openmc/lib/math.py @@ -199,7 +199,8 @@ def rotate_angle(uvw0, mu, phi, prn_seed=None): phi : float Azimuthal angle; if None, one will be sampled uniformly prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -231,7 +232,8 @@ def maxwell_spectrum(T, prn_seed=None): T : float Spectrum parameter prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -256,7 +258,8 @@ def watt_spectrum(a, b, prn_seed=None): b : float Spectrum parameter b prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns ------- @@ -281,7 +284,8 @@ def normal_variate(mean_value, std_dev, prn_seed=None): std_dev : float Standard deviation of the normal distribution prn_seed : int - PRNG seed; if None, one will be generated randomly + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. Returns -------