From 187160d0822ef9f4ceb23e489c96fd62336a1fb2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 19 Jan 2024 09:01:18 -0600 Subject: [PATCH] Ability to compute material volume fractions over mesh elements (#2802) --- include/openmc/capi.h | 3 + include/openmc/mesh.h | 30 +++++++- include/openmc/volume_calc.h | 36 +++++++++- openmc/lib/core.py | 5 +- openmc/lib/mesh.py | 84 +++++++++++++++++++++- src/mesh.cpp | 136 ++++++++++++++++++++++++++++++++++- src/volume_calc.cpp | 29 +------- tests/unit_tests/test_lib.py | 76 ++++++++++++++++++-- 8 files changed, 360 insertions(+), 39 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 1d4dc1a4ea9..063ed8d9fad 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -105,6 +105,9 @@ int openmc_mesh_filter_get_translation(int32_t index, double translation[3]); int openmc_mesh_filter_set_translation(int32_t index, double translation[3]); int openmc_mesh_get_id(int32_t index, int32_t* id); int openmc_mesh_set_id(int32_t index, int32_t id); +int openmc_mesh_get_n_elements(int32_t index, size_t* n); +int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin, + int result_size, void* result, int* hits, uint64_t* seed); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_new_filter(const char* type, int32_t* index); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index c88ffc3504c..37b77183248 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -9,6 +9,7 @@ #include "hdf5.h" #include "pugixml.hpp" #include "xtensor/xtensor.hpp" +#include #include "openmc/error.h" #include "openmc/memory.h" // for unique_ptr @@ -69,6 +70,12 @@ extern const libMesh::Parallel::Communicator* libmesh_comm; class Mesh { public: + // Types, aliases + struct MaterialVolume { + int32_t material; //!< material index + double volume; //!< volume in [cm^3] + }; + // Constructors and destructor Mesh() = default; Mesh(pugi::xml_node node); @@ -159,9 +166,28 @@ class Mesh { virtual std::string get_mesh_type() const = 0; + //! Determine volume of materials within a single mesh elemenet + // + //! \param[in] n_sample Number of samples within each element + //! \param[in] bin Index of mesh element + //! \param[out] Array of (material index, volume) for desired element + //! \param[inout] seed Pseudorandom number seed + //! \return Number of materials within element + int material_volumes(int n_sample, int bin, gsl::span volumes, + uint64_t* seed) const; + + //! Determine volume of materials within a single mesh elemenet + // + //! \param[in] n_sample Number of samples within each element + //! \param[in] bin Index of mesh element + //! \param[inout] seed Pseudorandom number seed + //! \return Vector of (material index, volume) for desired element + vector material_volumes( + int n_sample, int bin, uint64_t* seed) const; + // Data members - int id_ {-1}; //!< User-specified ID - int n_dimension_; //!< Number of dimensions + int id_ {-1}; //!< User-specified ID + int n_dimension_ {-1}; //!< Number of dimensions }; class StructuredMesh : public Mesh { diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 2773a39fa69..376b8c707dd 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -1,18 +1,23 @@ #ifndef OPENMC_VOLUME_CALC_H #define OPENMC_VOLUME_CALC_H +#include // for find #include #include +#include #include "openmc/array.h" +#include "openmc/openmp_interface.h" #include "openmc/position.h" #include "openmc/tallies/trigger.h" #include "openmc/vector.h" #include "pugixml.hpp" #include "xtensor/xtensor.hpp" - #include +#ifdef _OPENMP +#include +#endif namespace openmc { @@ -89,6 +94,35 @@ extern vector volume_calcs; // Non-member functions //============================================================================== +//! Reduce vector of indices and hits from each thread to a single copy +// +//! \param[in] local_indices Indices specific to each thread +//! \param[in] local_hits Hit count specific to each thread +//! \param[out] indices Reduced vector of indices +//! \param[out] hits Reduced vector of hits +template +void reduce_indices_hits(const vector& local_indices, + const vector& local_hits, vector& indices, vector& hits) +{ + const int n_threads = num_threads(); + +#pragma omp for ordered schedule(static, 1) + for (int i = 0; i < n_threads; ++i) { +#pragma omp ordered + for (int j = 0; j < local_indices.size(); ++j) { + // Check if this material has been added to the master list and if + // so, accumulate the number of hits + auto it = std::find(indices.begin(), indices.end(), local_indices[j]); + if (it == indices.end()) { + indices.push_back(local_indices[j]); + hits.push_back(local_hits[j]); + } else { + hits[it - indices.begin()] += local_hits[j]; + } + } + } +} + void free_memory_volume(); } // namespace openmc diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 25dc65fe84e..4ea1c86d0ba 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -611,7 +611,7 @@ def __set__(self, instance, value): class _FortranObject: def __repr__(self): - return "{}[{}]".format(type(self).__name__, self._index) + return "<{}(index={})>".format(type(self).__name__, self._index) class _FortranObjectWithID(_FortranObject): @@ -622,6 +622,9 @@ def __init__(self, uid=None, new=True, index=None): # OutOfBoundsError will be raised here by virtue of referencing self.id self.id + def __repr__(self): + return "<{}(id={})>".format(type(self).__name__, self.id) + @contextmanager def quiet_dll(output=True): diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 72d61cc6650..8a9953754ee 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -1,6 +1,8 @@ from collections.abc import Mapping -from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, - create_string_buffer) +from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure, + create_string_buffer, c_uint64, c_size_t) +from random import getrandbits +from typing import Optional, List, Tuple from weakref import WeakValueDictionary import numpy as np @@ -10,9 +12,18 @@ from . import _dll from .core import _FortranObjectWithID from .error import _error_handler +from .material import Material __all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'UnstructuredMesh', 'meshes'] + +class _MaterialVolume(Structure): + _fields_ = [ + ("material", c_int32), + ("volume", c_double) + ] + + # Mesh functions _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32), POINTER(c_int32)] @@ -24,6 +35,14 @@ _dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32] _dll.openmc_mesh_set_id.restype = c_int _dll.openmc_mesh_set_id.errcheck = _error_handler +_dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)] +_dll.openmc_mesh_get_n_elements.restype = c_int +_dll.openmc_mesh_get_n_elements.errcheck = _error_handler +_dll.openmc_mesh_material_volumes.argtypes = [ + c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume), + POINTER(c_int), POINTER(c_uint64)] +_dll.openmc_mesh_material_volumes.restype = c_int +_dll.openmc_mesh_material_volumes.errcheck = _error_handler _dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_mesh_index.restype = c_int _dll.openmc_get_mesh_index.errcheck = _error_handler @@ -123,6 +142,67 @@ def id(self): def id(self, mesh_id): _dll.openmc_mesh_set_id(self._index, mesh_id) + @property + def n_elements(self): + n = c_size_t() + _dll.openmc_mesh_get_n_elements(self._index, n) + return n.value + + def material_volumes( + self, + n_samples: int = 10_000, + prn_seed: Optional[int] = None + ) -> List[List[Tuple[Material, float]]]: + """Determine volume of materials in each mesh element + + .. versionadded:: 0.14.1 + + Parameters + ---------- + n_samples : int + Number of samples in each mesh element + prn_seed : int + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. + + Returns + ------- + List of tuple of (material, volume) for each mesh element. Void volume + is represented by having a value of None in the first element of a + tuple. + + """ + if n_samples <= 0: + raise ValueError("Number of samples must be positive") + if prn_seed is None: + prn_seed = getrandbits(63) + prn_seed = c_uint64(prn_seed) + + # Preallocate space for MaterialVolume results + size = 16 + result = (_MaterialVolume * size)() + + hits = c_int() # Number of materials hit in a given element + volumes = [] + for i_element in range(self.n_elements): + while True: + try: + _dll.openmc_mesh_material_volumes( + self._index, n_samples, i_element, size, result, hits, prn_seed) + except AllocationError: + # Increase size of result array and try again + size *= 2 + result = (_MaterialVolume * size)() + else: + # If no error, break out of loop + break + + volumes.append([ + (Material(index=r.material), r.volume) + for r in result[:hits.value] + ]) + return volumes + class RegularMesh(Mesh): """RegularMesh stored internally. diff --git a/src/mesh.cpp b/src/mesh.cpp index cb6851d04da..b70692347fe 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -22,15 +22,19 @@ #include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/geometry.h" #include "openmc/hdf5_interface.h" +#include "openmc/material.h" #include "openmc/memory.h" #include "openmc/message_passing.h" #include "openmc/openmp_interface.h" +#include "openmc/particle_data.h" #include "openmc/random_dist.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" +#include "openmc/volume_calc.h" #include "openmc/xml_interface.h" #ifdef LIBMESH @@ -141,6 +145,92 @@ vector Mesh::volumes() const return volumes; } +int Mesh::material_volumes( + int n_sample, int bin, gsl::span result, uint64_t* seed) const +{ + vector materials; + vector hits; + +#pragma omp parallel + { + vector local_materials; + vector local_hits; + GeometryState geom; + +#pragma omp for + for (int i = 0; i < n_sample; ++i) { + // Get seed for i-th sample + uint64_t seed_i = future_seed(3 * i, *seed); + + // Sample position and set geometry state + geom.r() = this->sample_element(bin, &seed_i); + geom.u() = {1., 0., 0.}; + geom.n_coord() = 1; + + // If this location is not in the geometry at all, move on to next block + if (!exhaustive_find_cell(geom)) + continue; + + int i_material = geom.material(); + + // Check if this material was previously hit and if so, increment count + auto it = + std::find(local_materials.begin(), local_materials.end(), i_material); + if (it == local_materials.end()) { + local_materials.push_back(i_material); + local_hits.push_back(1); + } else { + local_hits[it - local_materials.begin()]++; + } + } // omp for + + // Reduce index/hits lists from each thread into a single copy + reduce_indices_hits(local_materials, local_hits, materials, hits); + } // omp parallel + + // Advance RNG seed + advance_prn_seed(3 * n_sample, seed); + + // Make sure span passed in is large enough + if (hits.size() > result.size()) { + return -1; + } + + // Convert hits to fractions + for (int i_mat = 0; i_mat < hits.size(); ++i_mat) { + double fraction = double(hits[i_mat]) / n_sample; + result[i_mat].material = materials[i_mat]; + result[i_mat].volume = fraction * this->volume(bin); + } + return hits.size(); +} + +vector Mesh::material_volumes( + int n_sample, int bin, uint64_t* seed) const +{ + // Create result vector with space for 8 pairs + vector result; + result.reserve(8); + + int size = -1; + while (true) { + // Get material volumes + size = this->material_volumes( + n_sample, bin, {result.data(), result.data() + result.capacity()}, seed); + + // If capacity was sufficient, resize the vector and return + if (size >= 0) { + result.resize(size); + break; + } + + // Otherwise, increase capacity of the vector + result.reserve(2 * result.capacity()); + } + + return result; +} + //============================================================================== // Structured Mesh implementation //============================================================================== @@ -736,7 +826,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} fatal_error("Must specify either or on a mesh."); } - // Set volume fraction + // Set material volumes volume_frac_ = 1.0 / xt::prod(shape)(); element_volume_ = 1.0; @@ -1035,7 +1125,7 @@ double RectilinearMesh::volume(const MeshIndex& ijk) const double vol {1.0}; for (int i = 0; i < n_dimension_; i++) { - vol *= grid_[i][ijk[i] + 1] - grid_[i][ijk[i]]; + vol *= grid_[i][ijk[i]] - grid_[i][ijk[i] - 1]; } return vol; } @@ -1781,6 +1871,33 @@ extern "C" int openmc_mesh_set_id(int32_t index, int32_t id) return 0; } +//! Get the number of elements in a mesh +extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n) +{ + if (int err = check_mesh(index)) + return err; + *n = model::meshes[index]->n_bins(); + return 0; +} + +extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample, + int bin, int result_size, void* result, int* hits, uint64_t* seed) +{ + auto result_ = reinterpret_cast(result); + if (!result_) { + set_errmsg("Invalid result pointer passed to openmc_mesh_material_volumes"); + return OPENMC_E_INVALID_ARGUMENT; + } + + if (int err = check_mesh(index)) + return err; + + int n = model::meshes[index]->material_volumes( + n_sample, bin, {result_, result_ + result_size}, seed); + *hits = n; + return (n == -1) ? OPENMC_E_ALLOCATE : 0; +} + //! Get the dimension of a regular mesh extern "C" int openmc_regular_mesh_get_dimension( int32_t index, int** dims, int* n) @@ -1835,6 +1952,11 @@ extern "C" int openmc_regular_mesh_set_params( return err; RegularMesh* m = dynamic_cast(model::meshes[index].get()); + if (m->n_dimension_ == -1) { + set_errmsg("Need to set mesh dimension before setting parameters."); + return OPENMC_E_UNASSIGNED; + } + vector shape = {static_cast(n)}; if (ll && ur) { m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); @@ -1853,6 +1975,16 @@ extern "C" int openmc_regular_mesh_set_params( return OPENMC_E_INVALID_ARGUMENT; } + // Set material volumes + + // TODO: incorporate this into method in RegularMesh that can be called from + // here and from constructor + m->volume_frac_ = 1.0 / xt::prod(m->get_x_shape())(); + m->element_volume_ = 1.0; + for (int i = 0; i < m->n_dimension_; i++) { + m->element_volume_ *= m->width_[i]; + } + return 0; } diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 6a47e9490b5..8b5c27f14eb 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -203,32 +203,9 @@ vector VolumeCalculation::execute() const // At this point, each thread has its own pair of index/hits lists and we // now need to reduce them. OpenMP is not nearly smart enough to do this // on its own, so we have to manually reduce them - const int n_threads = num_threads(); - -#pragma omp for ordered schedule(static) - for (int i = 0; i < n_threads; ++i) { -#pragma omp ordered - for (int i_domain = 0; i_domain < n; ++i_domain) { - for (int j = 0; j < indices[i_domain].size(); ++j) { - // Check if this material has been added to the master list and if - // so, accumulate the number of hits - bool already_added = false; - for (int k = 0; k < master_indices[i_domain].size(); k++) { - if (indices[i_domain][j] == master_indices[i_domain][k]) { - master_hits[i_domain][k] += hits[i_domain][j]; - already_added = true; - break; - } - } - if (!already_added) { - // If we made it here, the material hasn't yet been added to the - // master list, so add entries to the master indices and master - // hits lists - master_indices[i_domain].push_back(indices[i_domain][j]); - master_hits[i_domain].push_back(hits[i_domain][j]); - } - } - } + for (int i_domain = 0; i_domain < n; ++i_domain) { + reduce_indices_hits(indices[i_domain], hits[i_domain], + master_indices[i_domain], master_hits[i_domain]); } } // omp parallel diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 11d6b12c9f9..8aa2ea6e9cd 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -1,4 +1,5 @@ from collections.abc import Mapping +from math import pi import os import numpy as np @@ -126,7 +127,7 @@ def test_cell(lib_init): cell = openmc.lib.cells[1] assert isinstance(cell.fill, openmc.lib.Material) cell.fill = openmc.lib.materials[1] - assert str(cell) == 'Cell[0]' + assert str(cell) == '' assert cell.name == "Fuel" cell.name = "Not fuel" assert cell.name == "Not fuel" @@ -584,6 +585,25 @@ def test_regular_mesh(lib_init): msf.translation = translation assert msf.translation == translation + # Test material volumes + mesh = openmc.lib.RegularMesh() + mesh.dimension = (2, 2, 1) + mesh.set_parameters(lower_left=(-0.63, -0.63, -0.5), + upper_right=(0.63, 0.63, 0.5)) + vols = mesh.material_volumes() + assert len(vols) == 4 + for elem_vols in vols: + assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26 / 4) + + # If the mesh extends beyond the boundaries of the model, the volumes should + # still be reported correctly + mesh.dimension = (1, 1, 1) + mesh.set_parameters(lower_left=(-1.0, -1.0, -0.5), + upper_right=(1.0, 1.0, 0.5)) + vols = mesh.material_volumes(100_000) + for elem_vols in vols: + assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26, 1e-2) + def test_rectilinear_mesh(lib_init): mesh = openmc.lib.RectilinearMesh() @@ -604,7 +624,7 @@ def test_rectilinear_mesh(lib_init): meshes = openmc.lib.meshes assert isinstance(meshes, Mapping) - assert len(meshes) == 2 + assert len(meshes) == 3 mesh = meshes[mesh.id] assert isinstance(mesh, openmc.lib.RectilinearMesh) @@ -615,8 +635,21 @@ def test_rectilinear_mesh(lib_init): msf = openmc.lib.MeshSurfaceFilter(mesh) assert msf.mesh == mesh + # Test material volumes + mesh = openmc.lib.RectilinearMesh() + w = 1.26 + mesh.set_grid([-w/2, -w/4, w/2], [-w/2, -w/4, w/2], [-0.5, 0.5]) + + vols = mesh.material_volumes() + assert len(vols) == 4 + assert sum(f[1] for f in vols[0]) == pytest.approx(w/4 * w/4) + assert sum(f[1] for f in vols[1]) == pytest.approx(w/4 * 3*w/4) + assert sum(f[1] for f in vols[2]) == pytest.approx(3*w/4 * w/4) + assert sum(f[1] for f in vols[3]) == pytest.approx(3*w/4 * 3*w/4) + + def test_cylindrical_mesh(lib_init): - deg2rad = lambda deg: deg*np.pi/180 + deg2rad = lambda deg: deg*pi/180 mesh = openmc.lib.CylindricalMesh() r_grid = [0., 5., 10.] phi_grid = np.radians([0., 10., 20.]) @@ -635,7 +668,7 @@ def test_cylindrical_mesh(lib_init): meshes = openmc.lib.meshes assert isinstance(meshes, Mapping) - assert len(meshes) == 3 + assert len(meshes) == 5 mesh = meshes[mesh.id] assert isinstance(mesh, openmc.lib.CylindricalMesh) @@ -646,6 +679,21 @@ def test_cylindrical_mesh(lib_init): msf = openmc.lib.MeshSurfaceFilter(mesh) assert msf.mesh == mesh + # Test material volumes + mesh = openmc.lib.CylindricalMesh() + r_grid = (0., 0.25, 0.5) + phi_grid = np.linspace(0., 2.0*pi, 4) + z_grid = (-0.5, 0.5) + mesh.set_grid(r_grid, phi_grid, z_grid) + + vols = mesh.material_volumes() + assert len(vols) == 6 + for i in range(0, 6, 2): + assert sum(f[1] for f in vols[i]) == pytest.approx(pi * 0.25**2 / 3) + for i in range(1, 6, 2): + assert sum(f[1] for f in vols[i]) == pytest.approx(pi * (0.5**2 - 0.25**2) / 3) + + def test_spherical_mesh(lib_init): deg2rad = lambda deg: deg*np.pi/180 mesh = openmc.lib.SphericalMesh() @@ -666,7 +714,7 @@ def test_spherical_mesh(lib_init): meshes = openmc.lib.meshes assert isinstance(meshes, Mapping) - assert len(meshes) == 4 + assert len(meshes) == 7 mesh = meshes[mesh.id] assert isinstance(mesh, openmc.lib.SphericalMesh) @@ -677,6 +725,24 @@ def test_spherical_mesh(lib_init): msf = openmc.lib.MeshSurfaceFilter(mesh) assert msf.mesh == mesh + # Test material volumes + mesh = openmc.lib.SphericalMesh() + r_grid = (0., 0.25, 0.5) + theta_grid = np.linspace(0., pi, 3) + phi_grid = np.linspace(0., 2.0*pi, 4) + mesh.set_grid(r_grid, theta_grid, phi_grid) + + vols = mesh.material_volumes() + assert len(vols) == 12 + d_theta = theta_grid[1] - theta_grid[0] + d_phi = phi_grid[1] - phi_grid[0] + for i in range(0, 12, 2): + assert sum(f[1] for f in vols[i]) == pytest.approx( + 0.25**3 / 3 * d_theta * d_phi * 2/pi) + for i in range(1, 12, 2): + assert sum(f[1] for f in vols[i]) == pytest.approx( + (0.5**3 - 0.25**3) / 3 * d_theta * d_phi * 2/pi) + def test_restart(lib_init, mpi_intracomm): # Finalize and re-init to make internal state consistent with XML.