diff --git a/Dockerfile b/Dockerfile index 341390c9a..d408b9de5 100644 --- a/Dockerfile +++ b/Dockerfile @@ -15,22 +15,24 @@ # sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number -FROM debian:bullseye-slim AS dependencies + +# global ARG as these ARGS are used in multiple stages +# By default one core is used to compile +ARG compile_cores=1 # By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support ARG build_dagmc=off ARG build_libmesh=off -# By default one core is used to compile -ARG compile_cores=1 +FROM debian:bullseye-slim AS dependencies + +ARG compile_cores +ARG build_dagmc +ARG build_libmesh # Set default value of HOME to /root ENV HOME=/root -# OpenMC variables -ARG openmc_branch=master -ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' - # Embree variables ENV EMBREE_TAG='v3.12.2' ENV EMBREE_REPO='https://github.com/embree/embree' @@ -60,7 +62,6 @@ ENV NJOY_REPO='https://github.com/njoy/NJOY2016' # Setup environment variables for Docker image ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \ - OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \ OPENMC_ENDF_DATA=/root/endf-b-vii.1 \ DEBIAN_FRONTEND=noninteractive @@ -174,12 +175,25 @@ RUN if [ "$build_libmesh" = "on" ]; then \ FROM dependencies AS build +ENV HOME=/root + +ARG openmc_branch=master +ENV OPENMC_REPO='https://github.com/openmc-dev/openmc' + +ARG compile_cores +ARG build_dagmc +ARG build_libmesh + +ENV DAGMC_INSTALL_DIR=$HOME/DAGMC/ +ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH + # clone and install openmc RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ && git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \ && mkdir build && cd build ; \ if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "on" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_DAGMC=on \ @@ -188,6 +202,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "off" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_DAGMC=ON \ @@ -195,6 +210,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "on" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on \ -DOPENMC_USE_LIBMESH=on \ @@ -202,6 +218,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ fi ; \ if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "off" ]; then \ cmake ../openmc \ + -DCMAKE_CXX_COMPILER=mpicxx \ -DOPENMC_USE_MPI=on \ -DHDF5_PREFER_PARALLEL=on ; \ fi ; \ @@ -211,5 +228,8 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \ FROM build AS release +ENV HOME=/root +ENV OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml + # Download cross sections (NNDC and WMP) and ENDF data needed by test suite RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index ef61604cb..8174108e1 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -339,6 +339,16 @@ Incoherent elastic scattering [eV\ :math:`^{-1}`]. :Attributes: - **type** (*char[]*) -- 'IncoherentElastic' +Sum of functions +---------------- + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "Sum" + - **n** (*int*) -- Number of functions +:Datasets: + - ***func_** (:ref:`function <1d_functions>`) -- Dataset for the + i-th function (indexing starts at 1) + .. _angle_energy: -------------------------- @@ -501,6 +511,19 @@ equiprobable bins. - **skewed** (*int8_t*) -- Whether discrete angles are equi-probable (0) or have a skewed distribution (1). +Mixed Elastic +------------- + +This angle-energy distribution is used when an evaluation specifies both +coherent and incoherent elastic thermal neutron scattering. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "mixed_elastic" +:Groups: - **coherent** -- Distribution for coherent elastic scattering. The + format is given in :ref:`angle_energy`. + - **incoherent** -- Distribution for incoherent elastic scattering. + The format is given in :ref:`angle_energy`. + .. _energy_distribution: -------------------- diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index f96d7181d..c1598f2c3 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -79,9 +79,15 @@ The current version of the statepoint file format is 17.0. - **width** (*double[]*) -- Width of each mesh cell in each dimension. - **Unstructured Mesh Only:** + - **filename** (*char[]*) -- Name of the mesh file. + - **library** (*char[]*) -- Mesh library used to represent the + mesh ("moab" or "libmesh"). + - **length_multiplier** (*double*) Scaling factor applied to the mesh. - **volumes** (*double[]*) -- Volume of each mesh cell. - - **centroids** (*double[]*) -- Location of the mesh cell - centroids. + - **vertices** (*double[]*) -- x, y, z values of the mesh vertices. + - **connectivity** (*int[]*) -- Connectivity array for the mesh + cells. + - **element_types** (*int[]*) -- Mesh element types. **/tallies/filters/** diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 03be25abf..882e3c71b 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/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 287738774..c74769be6 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -116,6 +116,7 @@ Angle-Energy Distributions IncoherentElasticAE IncoherentElasticAEDiscrete IncoherentInelasticAEDiscrete + MixedElasticAE Resonance Data -------------- diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 0929a11f7..69d3ff1f6 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/include/openmc/constants.h b/include/openmc/constants.h index 73f96ff95..a7362ea9e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -24,7 +24,7 @@ using double_4dvec = vector>>>; constexpr int HDF5_VERSION[] {3, 0}; // Version numbers for binary files -constexpr array VERSION_STATEPOINT {17, 0}; +constexpr array VERSION_STATEPOINT {18, 0}; constexpr array VERSION_PARTICLE_RESTART {2, 0}; constexpr array VERSION_TRACK {3, 0}; constexpr array VERSION_SUMMARY {6, 0}; diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 7beb8e452..e580874b4 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -126,6 +126,26 @@ private: debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] }; +//============================================================================== +//! Sum of multiple 1D functions +//============================================================================== + +class Sum1D : public Function1D { +public: + // Constructors + explicit Sum1D(hid_t group); + + //! Evaluate each function and sum results + //! \param[in] x independent variable + //! \return Function evaluated at x + double operator()(double E) const override; + + const unique_ptr& functions(int i) const { return functions_[i]; } + +private: + vector> functions_; //!< individual functions +}; + //! Read 1D function from HDF5 dataset //! \param[in] group HDF5 group containing dataset //! \param[in] name Name of dataset diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index e9d9b4f96..0092c08f8 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false); vector group_names(hid_t group_id); vector object_shape(hid_t obj_id); std::string object_name(hid_t obj_id); +hid_t open_object(hid_t group_id, const std::string& name); +void close_object(hid_t obj_id); //============================================================================== // Fortran compatibility functions diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 21c58862d..eece0f1ff 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -7,13 +7,14 @@ #include #include "hdf5.h" -#include "pugixml.hpp" #include "xtensor/xtensor.hpp" +#include "pugixml.hpp" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" #include "openmc/position.h" #include "openmc/vector.h" +#include "openmc/xml_interface.h" #ifdef DAGMC #include "moab/AdaptiveKDTree.hpp" @@ -36,6 +37,12 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + +enum class ElementType { UNSUPPORTED=-1, LINEAR_TET, LINEAR_HEX }; + //============================================================================== // Global variables //============================================================================== @@ -53,7 +60,7 @@ extern vector> meshes; #ifdef LIBMESH namespace settings { -// used when creating new libMesh::Mesh instances +// used when creating new libMesh::MeshBase instances extern unique_ptr libmesh_init; extern const libMesh::Parallel::Communicator* libmesh_comm; } // namespace settings @@ -485,6 +492,23 @@ public: //! \return The centroid of the bin virtual Position centroid(int bin) const = 0; + //! Get the number of vertices in the mesh + // + //! \return Number of vertices + virtual int n_vertices() const = 0; + + //! Retrieve a vertex of the mesh + // + //! \param[in] vertex ID + //! \return vertex coordinates + virtual Position vertex(int id) const = 0; + + //! Retrieve connectivity of a mesh element + // + //! \param[in] element ID + //! \return element connectivity as IDs of the vertices + virtual std::vector connectivity(int id) const = 0; + //! Get the volume of a mesh bin // //! \param[in] bin Bin to return the volume for @@ -557,6 +581,12 @@ public: Position centroid(int bin) const override; + int n_vertices() const override; + + Position vertex(int id) const override; + + std::vector connectivity(int id) const override; + double volume(int bin) const override; private: @@ -621,6 +651,9 @@ private: //! \return MOAB EntityHandle of tet moab::EntityHandle get_ent_handle_from_bin(int bin) const; + //! Get a vertex index into the global range from a handle + int get_vert_idx_from_handle(moab::EntityHandle vert) const; + //! Get the bin for a given mesh cell index // //! \param[in] idx Index of the mesh cell. @@ -655,6 +688,7 @@ private: // Data members moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh + moab::Range verts_; //!< Range of vertex EntityHandle's in the mesh moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree std::shared_ptr mbi_; //!< MOAB instance @@ -671,7 +705,8 @@ class LibMesh : public UnstructuredMesh { public: // Constructors LibMesh(pugi::xml_node node); - LibMesh(const std::string& filename, double length_multiplier = 1.0); + LibMesh(const std::string & filename, double length_multiplier = 1.0); + LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier = 1.0); static const std::string mesh_lib_type; @@ -701,10 +736,19 @@ public: Position centroid(int bin) const override; + int n_vertices() const override; + + Position vertex(int id) const override; + + std::vector connectivity(int id) const override; + double volume(int bin) const override; + libMesh::MeshBase* mesh_ptr() const { return m_; }; + private: void initialize() override; + void set_mesh_pointer_from_filename(const std::string& filename); // Methods @@ -715,7 +759,8 @@ private: int get_bin_from_element(const libMesh::Elem* elem) const; // Data members - unique_ptr m_; //!< pointer to the libMesh mesh instance + unique_ptr unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC + libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization vector> pl_; //!< per-thread point locators unique_ptr diff --git a/include/openmc/position.h b/include/openmc/position.h index 5ab0774f4..4fde420cc 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -5,6 +5,7 @@ #include #include // for out_of_range +#include "fmt/format.h" #include "openmc/array.h" #include "openmc/vector.h" @@ -210,4 +211,14 @@ using Direction = Position; } // namespace openmc +template<> +struct fmt::formatter : formatter { + template + auto format(const openmc::Position& pos, FormatContext& ctx) + { + return fmt::formatter::format( + fmt::format("({}, {}, {})", pos.x, pos.y, pos.z), ctx); + } +}; + #endif // OPENMC_POSITION_H diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 84ed68f51..5b18902af 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -150,6 +150,34 @@ private: //!< each incident energy }; +//============================================================================== +//! Mixed coherent/incoherent elastic angle-energy distribution +//============================================================================== + +class MixedElasticAE : public AngleEnergy { +public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group + explicit MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs); + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + //! \param[inout] seed Pseudorandom number seed pointer + void sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const override; + +private: + CoherentElasticAE coherent_dist_; //!< Coherent distribution + unique_ptr incoherent_dist_; //!< Incoherent distribution + + const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS + const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS +}; + } // namespace openmc #endif // OPENMC_SECONDARY_THERMAL_H diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 6009dc748..b8fde5478 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -44,6 +44,8 @@ class AngleEnergy(EqualityMixin, ABC): return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) elif dist_type == 'incoherent_inelastic': return openmc.data.IncoherentInelasticAE.from_hdf5(group) + elif dist_type == 'mixed_elastic': + return openmc.data.MixedElasticAE.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/function.py b/openmc/data/function.py index b5aa2117d..b0390d19c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -544,7 +544,7 @@ class Combination(EqualityMixin): self._operations = operations -class Sum(EqualityMixin): +class Sum(Function1D): """Sum of multiple functions. This class allows you to create a callable object which represents the sum @@ -578,6 +578,49 @@ class Sum(EqualityMixin): cv.check_type('functions', functions, Iterable, Callable) self._functions = functions + def to_hdf5(self, group, name='xy'): + """Write sum of functions to an HDF5 group + + .. versionadded:: 0.13.1 + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + sum_group = group.create_group(name) + sum_group.attrs['type'] = np.string_(type(self).__name__) + sum_group.attrs['n'] = len(self.functions) + for i, f in enumerate(self.functions): + f.to_hdf5(sum_group, f'func_{i+1}') + + @classmethod + def from_hdf5(cls, group): + """Generate sum of functions from an HDF5 group + + .. versionadded:: 0.13.1 + + Parameters + ---------- + group : h5py.Group + Group to read from + + Returns + ------- + openmc.data.Sum + Functions read from the group + + """ + n = group.attrs['n'] + functions = [ + Function1D.from_hdf5(group[f'func_{i+1}']) + for i in range(n) + ] + return cls(functions) + class Regions1D(EqualityMixin): r"""Piecewise composition of multiple functions. diff --git a/openmc/data/reconstruct.pyx b/openmc/data/reconstruct.pyx index f63a155b1..cd0bbc38b 100644 --- a/openmc/data/reconstruct.pyx +++ b/openmc/data/reconstruct.pyx @@ -299,8 +299,8 @@ def reconstruct_slbw(slbw, double E): # Determine shift and penetration at modified energy if slbw._competitive[i]: Ex = E + slbw.q_value[l]*(A + 1)/A - rhoc = slbw.channel_radius[l](Ex) - rhochat = slbw.scattering_radius[l](Ex) + rhoc = k*slbw.channel_radius[l](Ex) + rhochat = k*slbw.scattering_radius[l](Ex) P_c, S_c = penetration_shift(l, rhoc) if Ex < 0: P_c = 0 diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 4e0a30b64..3d5f2df64 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -19,12 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, isotopes from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .function import Tabulated1D, Function1D +from .function import Tabulated1D, Function1D, Sum from .njoy import make_ace_thermal from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, IncoherentElasticAEDiscrete, IncoherentInelasticAEDiscrete, - IncoherentInelasticAE) + IncoherentInelasticAE, MixedElasticAE) _THERMAL_NAMES = { @@ -694,29 +694,53 @@ class ThermalScattering(EqualityMixin): # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] - n_mu = ace.nxs[6] + 1 if idx != 0: - n_energy = int(ace.xss[idx]) - energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV - P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] - - if ace.nxs[5] == 4: + if ace.nxs[5] in (4, 5): # Coherent elastic - xs = CoherentElastic(energy, P*EV_PER_MEV) - distribution = CoherentElasticAE(xs) + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + coherent_xs = CoherentElastic(energy, P*EV_PER_MEV) + coherent_dist = CoherentElasticAE(coherent_xs) # Coherent elastic shouldn't have angular distributions listed + n_mu = ace.nxs[6] + 1 assert n_mu == 0 - else: - # Incoherent elastic - xs = Tabulated1D(energy, P) + + if ace.nxs[5] in (3, 5): + # Incoherent elastic scattering -- first determine if both + # incoherent and coherent are present (mixed) + mixed = (ace.nxs[5] == 5) + + # Get cross section values + idx = ace.jxs[7] if mixed else ace.jxs[4] + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + + incoherent_xs = Tabulated1D(energy, values) # Angular distribution + n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1 assert n_mu > 0 - idx = ace.jxs[6] + idx = ace.jxs[9] if mixed else ace.jxs[6] mu_out = ace.xss[idx:idx + n_energy * n_mu] mu_out.shape = (n_energy, n_mu) - distribution = IncoherentElasticAEDiscrete(mu_out) + incoherent_dist = IncoherentElasticAEDiscrete(mu_out) + + if ace.nxs[5] == 3: + xs = incoherent_xs + dist = incoherent_dist + elif ace.nxs[5] == 4: + xs = coherent_xs + dist = coherent_dist + else: + # Create mixed cross section -- note that coherent must come + # first due to assumption on C++ side + xs = Sum([coherent_xs, incoherent_xs]) + + # Create mixed distribution + distribution = MixedElasticAE(coherent_dist, incoherent_dist) table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution}) @@ -802,7 +826,7 @@ class ThermalScattering(EqualityMixin): # Replace ACE data with ENDF data rx, rx_endf = data.elastic, data_endf.elastic for t in temperatures: - if isinstance(rx_endf.xs[t], IncoherentElastic): + if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)): rx.xs[t] = rx_endf.xs[t] rx.distribution[t] = rx_endf.distribution[t] @@ -832,20 +856,14 @@ class ThermalScattering(EqualityMixin): # Read coherent/incoherent elastic data elastic = None if (7, 2) in ev.section: - xs = {} - distribution = {} - - file_obj = StringIO(ev.section[7, 2]) - lhtr = endf.get_head_record(file_obj)[2] - if lhtr == 1: - # coherent elastic - + # Define helper functions to avoid duplication + def get_coherent_elastic(file_obj): # Get structure factor at first temperature params, S = endf.get_tab1_record(file_obj) strT = _temperature_str(params[0]) n_temps = params[2] bragg_edges = S.x - xs[strT] = CoherentElastic(bragg_edges, S.y) + xs = {strT: CoherentElastic(bragg_edges, S.y)} distribution = {strT: CoherentElasticAE(xs[strT])} # Get structure factor for subsequent temperatures @@ -854,15 +872,34 @@ class ThermalScattering(EqualityMixin): strT = _temperature_str(params[0]) xs[strT] = CoherentElastic(bragg_edges, S) distribution[strT] = CoherentElasticAE(xs[strT]) + return xs, distribution - elif lhtr == 2: - # incoherent elastic + def get_incoherent_elastic(file_obj): params, W = endf.get_tab1_record(file_obj) bound_xs = params[0] + xs = {} + distribution = {} for T, debye_waller in zip(W.x, W.y): strT = _temperature_str(T) xs[strT] = IncoherentElastic(bound_xs, debye_waller) distribution[strT] = IncoherentElasticAE(debye_waller) + return xs, distribution + + file_obj = StringIO(ev.section[7, 2]) + lhtr = endf.get_head_record(file_obj)[2] + if lhtr == 1: + # coherent elastic + xs, distribution = get_coherent_elastic(file_obj) + elif lhtr == 2: + # incoherent elastic + xs, distribution = get_incoherent_elastic(file_obj) + elif lhtr == 3: + # mixed coherent / incoherent elastic + xs_c, dist_c = get_coherent_elastic(file_obj) + xs_i, dist_i = get_incoherent_elastic(file_obj) + assert sorted(xs_c) == sorted(xs_i) + xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c} + distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c} elastic = ThermalScatteringReaction(xs, distribution) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 05393e7bb..17a560092 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -2,6 +2,7 @@ import numpy as np from .angle_energy import AngleEnergy from .correlated import CorrelatedAngleEnergy +import openmc.data class CoherentElasticAE(AngleEnergy): @@ -43,7 +44,27 @@ class CoherentElasticAE(AngleEnergy): """ group.attrs['type'] = np.string_('coherent_elastic') - group['coherent_xs'] = group.parent['xs'] + self.coherent_xs.to_hdf5(group, 'coherent_xs') + + @classmethod + def from_hdf5(cls, group): + """Generate coherent elastic distribution from HDF5 data + + .. versionadded:: 0.13.1 + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.CoherentElasticAE + Coherent elastic distribution + + """ + coherent_xs = openmc.data.CoherentElastic.from_hdf5(group['coherent_xs']) + return cls(coherent_xs) class IncoherentElasticAE(AngleEnergy): @@ -101,7 +122,7 @@ class IncoherentElasticAE(AngleEnergy): Incoherent elastic distribution """ - return cls(group['debye_waller']) + return cls(group['debye_waller'][()]) class IncoherentElasticAEDiscrete(AngleEnergy): @@ -210,3 +231,62 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): class IncoherentInelasticAE(CorrelatedAngleEnergy): _name = 'incoherent_inelastic' + + +class MixedElasticAE(AngleEnergy): + """Secondary distribution for mixed coherent/incoherent thermal elastic + + .. versionadded:: 0.13.1 + + Parameters + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + Attributes + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + """ + def __init__(self, coherent, incoherent): + self.coherent = coherent + self.incoherent = incoherent + + def to_hdf5(self, group): + """Write mixed elastic distribution to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + group.attrs['type'] = np.string_('mixed_elastic') + coherent_group = group.create_group('coherent') + self.coherent.to_hdf5(coherent_group) + incoherent_group = group.create_group('incoherent') + self.incoherent.to_hdf5(incoherent_group) + + @classmethod + def from_hdf5(cls, group): + """Generate mixed thermal elastic distribution from HDF5 data + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.MixedElasticAE + Mixed thermal elastic distribution + + """ + coherent = AngleEnergy.from_hdf5(group['coherent']) + incoherent = AngleEnergy.from_hdf5(group['incoherent']) + return cls(coherent, incoherent) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index de5f4adf2..479e414ae 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, c_size_t) 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 = [c_size_t, 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,45 @@ def run(output=True): _dll.openmc_run() +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 + Pseudorandom number generator (PRNG) seed; if None, one will be + generated randomly. + + Returns + ------- + 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 + sites_array = (_SourceSite * n_samples)() + _dll.openmc_sample_external_source(c_size_t(n_samples), c_uint64(prn_seed), sites_array) + + # 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(): """Initialize simulation""" _dll.openmc_simulation_init() diff --git a/openmc/lib/math.py b/openmc/lib/math.py index bdb09b2d1..0f61450bf 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] @@ -200,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 ------- @@ -232,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 ------- @@ -240,10 +241,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)) @@ -257,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 ------- @@ -265,7 +267,7 @@ def watt_spectrum(a, b, prn_seed=None): Sampled outgoing energy """ - + if prn_seed is None: prn_seed = getrandbits(63) @@ -282,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 ------- @@ -290,7 +293,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/openmc/material.py b/openmc/material.py index a837b1f17..6c2641680 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -34,7 +34,9 @@ class Material(IDManagerMixin): To create a material, one should create an instance of this class, add nuclides or elements with :meth:`Material.add_nuclide` or :meth:`Material.add_element`, respectively, and set the total material - density with :meth:`Material.set_density()`. The material can then be + density with :meth:`Material.set_density()`. Alternatively, you can + use :meth:`Material.add_components()` to pass a dictionary + containing all the component information. The material can then be assigned to a cell using the :attr:`Cell.fill` attribute. Parameters @@ -403,6 +405,54 @@ class Material(IDManagerMixin): self._nuclides.append(NuclideTuple(nuclide, percent, percent_type)) + def add_components(self, components: dict, percent_type: str = 'ao'): + """ Add multiple elements or nuclides to a material + + .. versionadded:: 0.13.1 + + Parameters + ---------- + components : dict of str to float or dict + Dictionary mapping element or nuclide names to their atom or weight + percent. To specify enrichment of an element, the entry of + ``components`` for that element must instead be a dictionary + containing the keyword arguments as well as a value for + ``'percent'`` + percent_type : {'ao', 'wo'} + 'ao' for atom percent and 'wo' for weight percent + + Examples + -------- + >>> mat = openmc.Material() + >>> components = {'Li': {'percent': 1.0, + >>> 'enrichment': 60.0, + >>> 'enrichment_target': 'Li7'}, + >>> 'Fl': 1.0, + >>> 'Be6': 0.5} + >>> mat.add_components(components) + + """ + + for component, params in components.items(): + cv.check_type('component', component, str) + if isinstance(params, float): + params = {'percent': params} + + else: + cv.check_type('params', params, dict) + if 'percent' not in params: + raise ValueError("An entry in the dictionary does not have " + "a required key: 'percent'") + + params['percent_type'] = percent_type + + ## check if nuclide + if str.isdigit(component[-1]): + self.add_nuclide(component, **params) + else: # is element + kwargs = params + self.add_element(component, **params) + def remove_nuclide(self, nuclide: str): """Remove a nuclide from the material diff --git a/openmc/mesh.py b/openmc/mesh.py index 1fe5429d0..daff4dedc 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2,6 +2,7 @@ from abc import ABC, abstractmethod from collections.abc import Iterable from math import pi from numbers import Real, Integral +from pathlib import Path import warnings from xml.etree import ElementTree as ET @@ -248,8 +249,14 @@ class StructuredMesh(MeshBase): # create VTK arrays for each of # the data sets + + # maintain a list of the datasets as added + # to the VTK arrays to ensure they persist + # in memory until the file is written + datasets_out = [] for label, dataset in datasets.items(): dataset = np.asarray(dataset).flatten() + datasets_out.append(dataset) if volume_normalization: dataset /= self.volumes.flatten() @@ -692,7 +699,7 @@ class RegularMesh(StructuredMesh): root_cell.fill = lattice return root_cell, cells - + def write_data_to_vtk(self, filename, datasets, volume_normalization=True): """Creates a VTK object of the mesh @@ -1434,7 +1441,7 @@ class UnstructuredMesh(MeshBase): Parameters ---------- - filename : str + filename : str or pathlib.Path Location of the unstructured mesh file library : {'moab', 'libmesh'} Mesh library used for the unstructured mesh tally @@ -1462,20 +1469,39 @@ class UnstructuredMesh(MeshBase): be generated for this mesh volumes : Iterable of float Volumes of the unstructured mesh elements + centroids : numpy.ndarray + Centroids of the mesh elements with array shape (n_elements, 3) + + vertices : numpy.ndarray + Coordinates of the mesh vertices with array shape (n_elements, 3) + + .. versionadded:: 0.13.1 + connectivity : numpy.ndarray + Connectivity of the elements with array shape (n_elements, 8) + + .. versionadded:: 0.13.1 + element_types : Iterable of integers + Mesh element types + + .. versionadded:: 0.13.1 total_volume : float Volume of the unstructured mesh in total - centroids : Iterable of tuple - An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0), - (1.0, 1.0, 1.0), ...] """ + + _UNSUPPORTED_ELEM = -1 + _LINEAR_TET = 0 + _LINEAR_HEX = 1 + def __init__(self, filename, library, mesh_id=None, name='', length_multiplier=1.0): super().__init__(mesh_id, name) self.filename = filename self._volumes = None - self._centroids = None + self._n_elements = None + self._conectivity = None + self._vertices = None self.library = library - self._output = True + self._output = False self.length_multiplier = length_multiplier @property @@ -1484,7 +1510,7 @@ class UnstructuredMesh(MeshBase): @filename.setter def filename(self, filename): - cv.check_type('Unstructured Mesh filename', filename, str) + cv.check_type('Unstructured Mesh filename', filename, (str, Path)) self._filename = filename @property @@ -1540,19 +1566,33 @@ class UnstructuredMesh(MeshBase): def centroids(self): return self._centroids + @property + def vertices(self): + return self._vertices + + @property + def connectivity(self): + return self._connectivity + + @property + def element_types(self): + return self._element_types + + @property + def centroids(self): + return np.array([self.centroid(i) for i in range(self.n_elements)]) + @property def n_elements(self): - if self._centroids is None: + if self._n_elements is None: raise RuntimeError("No information about this mesh has " "been loaded from a statepoint file.") - return len(self._centroids) + return self._n_elements - - @centroids.setter - def centroids(self, centroids): - cv.check_type("Unstructured mesh centroids", centroids, - Iterable, Real) - self._centroids = centroids + @n_elements.setter + def n_elements(self, val): + cv.check_type('Number of elements', val, Integral) + self._n_elements = val @property def length_multiplier(self): @@ -1567,29 +1607,44 @@ class UnstructuredMesh(MeshBase): @property def dimension(self): - return self.n_elements + return (self.n_elements,) @property def n_dimension(self): return 3 - @property - def vertices(self): - raise NotImplementedError("Vertices for UnstructuredMesh objects are " - "not yet available") - def __repr__(self): string = super().__repr__() string += '{: <16}=\t{}\n'.format('\tFilename', self.filename) - string += '{: <16}=\t{}\n'.format('\tMesh Library', self.mesh_lib) + string += '{: <16}=\t{}\n'.format('\tMesh Library', self.library) if self.length_multiplier != 1.0: string += '{: <16}=\t{}\n'.format('\tLength multiplier', self.length_multiplier) return string - def write_data_to_vtk(self, filename, datasets, volume_normalization=True): - """Map data to the unstructured mesh element centroids - to create a VTK point-cloud dataset. + def centroid(self, bin): + """Return the vertex averaged centroid of an element + + Parameters + ---------- + bin : int + Bin ID for the returned centroid + + Returns + ------- + numpy.ndarray + x, y, z values of the element centroid + + """ + conn = self.connectivity[bin] + coords = self.vertices[conn] + return coords.mean(axis=0) + + def write_vtk_mesh(self, **kwargs): + """Map data to unstructured VTK mesh elements. + + .. deprecated:: 0.13 + Use :func:`UnstructuredMesh.write_data_to_vtk` instead. Parameters ---------- @@ -1601,83 +1656,103 @@ class UnstructuredMesh(MeshBase): volume_normalization : bool Whether or not to normalize the data by the volume of the mesh elements - - Raises - ------ - RuntimeError - when the size of a dataset doesn't match the number of cells """ + warnings.warn( + "The 'UnstructuredMesh.write_vtk_mesh' method has been renamed " + "to 'write_data_to_vtk' and will be removed in a future version " + " of OpenMC.", FutureWarning + ) + self.write_data_to_vtk(**kwargs) + def write_data_to_vtk(self, filename=None, datasets=None, volume_normalization=True): + """Map data to unstructured VTK mesh elements. + + Parameters + ---------- + filename : str or pathlib.Path + Name of the VTK file to write + datasets : dict + Dictionary whose keys are the data labels + and values are numpy appropriately sized arrays + of the data + volume_normalization : bool + Whether or not to normalize the data by the + volume of the mesh elements + """ import vtk - from vtk.util import numpy_support as vtk_npsup + from vtk.util import numpy_support as nps - if self.centroids is None: - raise RuntimeError("No centroid information is present on this " - "unstructured mesh. Please load this " - "information from a relevant statepoint file.") + if self.connectivity is None or self.vertices is None: + raise RuntimeError('This mesh has not been ' + 'loaded from a statepoint file.') - if self.volumes is None and volume_normalization: - raise RuntimeError("No volume data is present on this " - "unstructured mesh. Please load the " - " mesh information from a statepoint file.") + if filename is None: + filename = f'mesh_{self.id}.vtk' - # check that the data sets are appropriately sized - errmsg = "The size of the dataset {} should be equal to the number of cells" - for label, dataset in datasets.items(): - if isinstance(dataset, np.ndarray): - if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) + writer = vtk.vtkUnstructuredGridWriter() + + writer.SetFileName(str(filename)) + + grid = vtk.vtkUnstructuredGrid() + + vtk_pnts = vtk.vtkPoints() + vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices)) + grid.SetPoints(vtk_pnts) + + n_skipped = 0 + elems = [] + for elem_type, conn in zip(self.element_types, self.connectivity): + if elem_type == self._LINEAR_TET: + elem = vtk.vtkTetra() + elif elem_type == self._LINEAR_HEX: + elem = vtk.vtkHexahedron() + elif elem_type == self._UNSUPPORTED_ELEM: + n_skipped += 1 else: - if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]: - raise RuntimeError(errmsg.format(label)) - cv.check_type('label', label, str) + raise RuntimeError(f'Invalid element type {elem_type} found') + for i, c in enumerate(conn): + if c == -1: + break + elem.GetPointIds().SetId(i, c) + elems.append(elem) - # create data arrays for the cells/points - cell_dim = 1 - vertices = vtk.vtkCellArray() - points = vtk.vtkPoints() + if n_skipped > 0: + warnings.warn(f'{n_skipped} elements were not written because ' + 'they are not of type linear tet/hex') - for centroid in self.centroids: - # create a point for each centroid - point_id = points.InsertNextPoint(centroid * self.length_multiplier) - # create a cell of type "Vertex" for each point - cell_id = vertices.InsertNextCell(cell_dim, (point_id,)) + for elem in elems: + grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) - # create a VTK data object - poly_data = vtk.vtkPolyData() - poly_data.SetPoints(points) - poly_data.SetVerts(vertices) - - # strange VTK nuance: - # data must be held in some container - # until the vtk file is written - data_holder = [] - - # create VTK arrays for each of - # the data sets - for label, dataset in datasets.items(): - dataset = np.asarray(dataset).flatten() + # check that datasets are the correct size + datasets_out = [] + if datasets is not None: + for name, data in datasets.items(): + if data.shape != self.dimension: + raise ValueError(f'Cannot apply dataset "{name}" with ' + f'shape {data.shape} to mesh {self.id} ' + f'with dimensions {self.dimension}') if volume_normalization: - dataset /= self.volumes.flatten() + for name, data in datasets.items(): + if np.issubdtype(data.dtype, np.integer): + warnings.warn(f'Integer data set "{name}" will ' + 'not be volume-normalized.') + continue + data /= self.volumes - array = vtk.vtkDoubleArray() - array.SetName(label) - array.SetNumberOfComponents(1) - array.SetArray(vtk_npsup.numpy_to_vtk(dataset), - dataset.size, - True) + # add data to the mesh + for name, data in datasets.items(): + datasets_out.append(data) + arr = vtk.vtkDoubleArray() + arr.SetName(name) + arr.SetNumberOfTuples(data.size) - data_holder.append(dataset) - poly_data.GetPointData().AddArray(array) + for i in range(data.size): + arr.SetTuple1(i, data.flat[i]) + grid.GetCellData().AddArray(arr) - # set filename - if not filename.endswith(".vtk"): - filename += ".vtk" + writer.SetInputData(grid) - writer = vtk.vtkGenericDataObjectWriter() - writer.SetFileName(filename) - writer.SetInputData(poly_data) writer.Write() @classmethod @@ -1688,10 +1763,14 @@ class UnstructuredMesh(MeshBase): mesh = cls(filename, library, mesh_id=mesh_id) vol_data = group['volumes'][()] - centroids = group['centroids'][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) - mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3)) - mesh.size = mesh.volumes.size + mesh.n_elements = mesh.volumes.size + + vertices = group['vertices'][()] + mesh._vertices = vertices.reshape((-1, 3)) + connectivity = group['connectivity'][()] + mesh._connectivity = connectivity.reshape((-1, 8)) + mesh._element_types = group['element_types'][()] if 'length_multiplier' in group: mesh.length_multiplier = group['length_multiplier'][()] @@ -1713,7 +1792,7 @@ class UnstructuredMesh(MeshBase): element.set("type", "unstructured") element.set("library", self._library) subelement = ET.SubElement(element, "filename") - subelement.text = self.filename + subelement.text = str(self.filename) if self._length_multiplier != 1.0: element.set("length_multiplier", str(self.length_multiplier)) diff --git a/openmc/mixin.py b/openmc/mixin.py index 516162464..31c26ec76 100644 --- a/openmc/mixin.py +++ b/openmc/mixin.py @@ -14,8 +14,11 @@ class EqualityMixin: def __eq__(self, other): if isinstance(other, type(self)): for key, value in self.__dict__.items(): - if not np.array_equal(value, other.__dict__.get(key)): - return False + if isinstance(value, np.ndarray): + if not np.array_equal(value, other.__dict__.get(key)): + return False + else: + return value == other.__dict__.get(key) else: return False diff --git a/openmc/region.py b/openmc/region.py index 56926187d..4e74a08ad 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -259,13 +259,16 @@ class Region(ABC): clone[:] = [n.clone(memo) for n in self] return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): """Translate region in given direction Parameters ---------- vector : iterable of float Direction in which region should be translated + inplace : bool + Whether or not to return a region based on new surfaces or one based + on the original surfaces that have been modified. memo : dict or None Dictionary used for memoization. This parameter is used internally and should not be specified by the user. @@ -279,7 +282,7 @@ class Region(ABC): if memo is None: memo = {} - return type(self)(n.translate(vector, memo) for n in self) + return type(self)(n.translate(vector, inplace, memo) for n in self) def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): @@ -308,7 +311,7 @@ class Region(ABC): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. memo : dict or None @@ -622,10 +625,10 @@ class Complement(Region): clone.node = self.node.clone(memo) return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): if memo is None: memo = {} - return type(self)(self.node.translate(vector, memo)) + return type(self)(self.node.translate(vector, inplace, memo)) def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): diff --git a/openmc/settings.py b/openmc/settings.py index 2d8121659..47c61e6ae 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -75,7 +75,7 @@ class Settings: Maximum number of lost particles .. versionadded:: 0.12 - rel_max_lost_particles : int + rel_max_lost_particles : float Maximum number of lost particles, relative to the total number of particles .. versionadded:: 0.12 @@ -298,191 +298,191 @@ class Settings: self._max_tracks = None @property - def run_mode(self): + def run_mode(self) -> str: return self._run_mode.value @property - def batches(self): + def batches(self) -> int: return self._batches @property - def generations_per_batch(self): + def generations_per_batch(self) -> int: return self._generations_per_batch @property - def inactive(self): + def inactive(self) -> int: return self._inactive @property - def max_lost_particles(self): + def max_lost_particles(self) -> int: return self._max_lost_particles @property - def rel_max_lost_particles(self): + def rel_max_lost_particles(self) -> float: return self._rel_max_lost_particles @property - def particles(self): + def particles(self) -> int: return self._particles @property - def keff_trigger(self): + def keff_trigger(self) -> dict: return self._keff_trigger @property - def energy_mode(self): + def energy_mode(self) -> str: return self._energy_mode @property - def max_order(self): + def max_order(self) -> int: return self._max_order @property - def source(self): + def source(self) -> typing.List[Source]: return self._source @property - def confidence_intervals(self): + def confidence_intervals(self) -> bool: return self._confidence_intervals @property - def electron_treatment(self): + def electron_treatment(self) -> str: return self._electron_treatment @property - def ptables(self): + def ptables(self) -> bool: return self._ptables @property - def photon_transport(self): + def photon_transport(self) -> bool: return self._photon_transport @property - def seed(self): + def seed(self) -> int: return self._seed @property - def survival_biasing(self): + def survival_biasing(self) -> bool: return self._survival_biasing @property - def entropy_mesh(self): + def entropy_mesh(self) -> RegularMesh: return self._entropy_mesh @property - def trigger_active(self): + def trigger_active(self) -> bool: return self._trigger_active @property - def trigger_max_batches(self): + def trigger_max_batches(self) -> int: return self._trigger_max_batches @property - def trigger_batch_interval(self): + def trigger_batch_interval(self) -> int: return self._trigger_batch_interval @property - def output(self): + def output(self) -> dict: return self._output @property - def sourcepoint(self): + def sourcepoint(self) -> dict: return self._sourcepoint @property - def statepoint(self): + def statepoint(self) -> dict: return self._statepoint @property - def surf_source_read(self): + def surf_source_read(self) -> dict: return self._surf_source_read @property - def surf_source_write(self): + def surf_source_write(self) -> dict: return self._surf_source_write @property - def no_reduce(self): + def no_reduce(self) -> bool: return self._no_reduce @property - def verbosity(self): + def verbosity(self) -> int: return self._verbosity @property - def tabular_legendre(self): + def tabular_legendre(self) -> dict: return self._tabular_legendre @property - def temperature(self): + def temperature(self) -> dict: return self._temperature @property - def trace(self): + def trace(self) -> typing.Iterable: return self._trace @property - def track(self): + def track(self) -> typing.Iterable[typing.Iterable[int]]: return self._track @property - def cutoff(self): + def cutoff(self) -> dict: return self._cutoff @property - def ufs_mesh(self): + def ufs_mesh(self) -> RegularMesh: return self._ufs_mesh @property - def resonance_scattering(self): + def resonance_scattering(self) -> dict: return self._resonance_scattering @property - def volume_calculations(self): + def volume_calculations(self) -> typing.List[VolumeCalculation]: return self._volume_calculations @property - def create_fission_neutrons(self): + def create_fission_neutrons(self) -> bool: return self._create_fission_neutrons @property - def delayed_photon_scaling(self): + def delayed_photon_scaling(self) -> bool: return self._delayed_photon_scaling @property - def material_cell_offsets(self): + def material_cell_offsets(self) -> bool: return self._material_cell_offsets @property - def log_grid_bins(self): + def log_grid_bins(self) -> int: return self._log_grid_bins @property - def event_based(self): + def event_based(self) -> bool: return self._event_based @property - def max_particles_in_flight(self): + def max_particles_in_flight(self) -> int: return self._max_particles_in_flight @property - def write_initial_source(self): + def write_initial_source(self) -> bool: return self._write_initial_source @property - def weight_windows(self): + def weight_windows(self) -> typing.List[WeightWindows]: return self._weight_windows @property - def weight_windows_on(self): + def weight_windows_on(self) -> bool: return self._weight_windows_on @property - def max_splits(self): + def max_splits(self) -> int: return self._max_splits @property - def max_tracks(self): + def max_tracks(self) -> int: return self._max_tracks @run_mode.setter @@ -517,7 +517,7 @@ class Settings: self._max_lost_particles = max_lost_particles @rel_max_lost_particles.setter - def rel_max_lost_particles(self, rel_max_lost_particles: int): + def rel_max_lost_particles(self, rel_max_lost_particles: float): cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real) cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0) cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 630675816..138c0c718 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -11,7 +11,7 @@ from uncertainties import ufloat import openmc import openmc.checkvalue as cv -_VERSION_STATEPOINT = 17 +_VERSION_STATEPOINT = 18 class StatePoint: diff --git a/openmc/surface.py b/openmc/surface.py index b3bcebb11..2996897a4 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -336,9 +336,9 @@ class Surface(IDManagerMixin, ABC): ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether or not to return a new instance of this Surface or to - modify the coefficients of this Surface. Defaults to False + modify the coefficients of this Surface. Returns ------- @@ -374,7 +374,7 @@ class Surface(IDManagerMixin, ABC): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. @@ -568,9 +568,9 @@ class PlaneMixin: ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether or not to return a new instance of a Plane or to modify the - coefficients of this plane. Defaults to False + coefficients of this plane. Returns ------- @@ -1012,9 +1012,8 @@ class QuadricMixin: ---------- vector : iterable of float Direction in which surface should be translated - inplace : boolean + inplace : bool Whether to return a clone of the Surface or the Surface itself. - Defaults to False Returns ------- @@ -2563,7 +2562,7 @@ class Halfspace(Region): clone.surface = self.surface.clone(memo) return clone - def translate(self, vector, memo=None): + def translate(self, vector, inplace=False, memo=None): """Translate half-space in given direction Parameters @@ -2585,7 +2584,7 @@ class Halfspace(Region): # If translated surface not in memo, add it key = (self.surface, tuple(vector)) if key not in memo: - memo[key] = self.surface.translate(vector) + memo[key] = self.surface.translate(vector, inplace) # Return translated half-space return type(self)(memo[key], self.side) @@ -2617,7 +2616,7 @@ class Halfspace(Region): :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. - inplace : boolean + inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. memo : dict or None diff --git a/openmc/universe.py b/openmc/universe.py index b65263ef9..78f75bbfc 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -12,9 +12,12 @@ import numpy as np import openmc import openmc.checkvalue as cv + from ._xml import get_text +from .checkvalue import check_type, check_value from .mixin import IDManagerMixin from .plots import _SVG_COLORS +from .surface import _BOUNDARY_TYPES class UniverseBase(ABC, IDManagerMixin): @@ -713,6 +716,78 @@ class DAGMCUniverse(UniverseBase): dagmc_element.set('filename', self.filename) xml_element.append(dagmc_element) + def bounding_region(self, bounded_type='box', boundary_type='vacuum'): + """Creates a either a spherical or box shaped bounding region around + the DAGMC geometry. + Parameters + ---------- + bounded_type : str + The type of bounding surface(s) to use when constructing the region. + Options include a single spherical surface (sphere) or a rectangle + made from six planes (box). + boundary_type : str + Boundary condition that defines the behavior for particles hitting + the surface. Defaults to vacuum boundary condition. Passed into the + surface construction. + Returns + ------- + openmc.Region + Region instance + """ + + check_type('boundary type', boundary_type, str) + check_value('boundary type', boundary_type, _BOUNDARY_TYPES) + check_type('bounded type', bounded_type, str) + check_value('bounded type', bounded_type, ('box', 'sphere')) + + bounding_box = self.bounding_box + + if bounded_type == 'sphere': + import math + bounding_box_center = (bounding_box[0] + bounding_box[1])/2 + radius = math.dist(bounding_box[0], bounding_box[1]) + bounding_surface = openmc.Sphere( + x0=bounding_box_center[0], + y0=bounding_box_center[1], + z0=bounding_box_center[2], + boundary_type=boundary_type, + r=radius, + ) + + return -bounding_surface + + if bounded_type == 'box': + # defines plane surfaces for all six faces of the bounding box + lower_x = openmc.XPlane(bounding_box[0][0], boundary_type=boundary_type) + upper_x = openmc.XPlane(bounding_box[1][0], boundary_type=boundary_type) + lower_y = openmc.YPlane(bounding_box[0][1], boundary_type=boundary_type) + upper_y = openmc.YPlane(bounding_box[1][1], boundary_type=boundary_type) + lower_z = openmc.ZPlane(bounding_box[0][2], boundary_type=boundary_type) + upper_z = openmc.ZPlane(bounding_box[1][2], boundary_type=boundary_type) + + return +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z + + def bounded_universe(self, bounding_cell_id=10000, **kwargs): + """Returns an openmc.Universe filled with this DAGMCUniverse and bounded + with a cell. Defaults to a box cell with a vacuum surface however this + can be changed using the kwargs which are passed directly to + DAGMCUniverse.bounding_region(). + + Parameters + ---------- + bounding_cell_id : int + The cell ID number to use for the bounding cell, defaults to 1000 to reduce + the chance of overlapping ID numbers with the DAGMC geometry. + + Returns + ------- + openmc.Universe + Universe instance + """ + + bounding_cell = openmc.Cell(fill=self, cell_id =bounding_cell_id, region=self.bounding_region(**kwargs)) + return openmc.Universe(cells=[bounding_cell]) + @classmethod def from_hdf5(cls, group): """Create DAGMC universe from HDF5 group diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 2e5b7d8e0..c5b7a00b6 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -291,10 +291,10 @@ class WeightWindows(IDManagerMixin): subelement.text = ' '.join(str(e) for e in self.energy_bounds) subelement = ET.SubElement(element, 'lower_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds) + subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel('F')) subelement = ET.SubElement(element, 'upper_ww_bounds') - subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds) + subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel('F')) subelement = ET.SubElement(element, 'survival_ratio') subelement.text = str(self.survival_ratio) diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 8ce5aabda..dc3bc2b2c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -183,7 +183,9 @@ void DAGUniverse::init_geometry() } else { warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3))); fatal_error(fmt::format("Cell ID {} exists in both DAGMC Universe {} " - "and the CSG geometry.", + "and the CSG geometry. Setting auto_geom_ids " + "to True when initiating the DAGMC Universe may " + "resolve this issue", c->id_, this->id_)); } diff --git a/src/endf.cpp b/src/endf.cpp index b42c8641d..c0c1d2e7e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -90,23 +90,25 @@ bool is_inelastic_scatter(int mt) unique_ptr read_function(hid_t group, const char* name) { - hid_t dset = open_dataset(group, name); + hid_t obj_id = open_object(group, name); std::string func_type; - read_attribute(dset, "type", func_type); + read_attribute(obj_id, "type", func_type); unique_ptr func; if (func_type == "Tabulated1D") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "Polynomial") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "CoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "IncoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); + } else if (func_type == "Sum") { + func = make_unique(obj_id); } else { throw std::runtime_error {"Unknown function type " + func_type + - " for dataset " + object_name(dset)}; + " for dataset " + object_name(obj_id)}; } - close_dataset(dset); + close_object(obj_id); return func; } @@ -271,4 +273,30 @@ double IncoherentElasticXS::operator()(double E) const return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); } +//============================================================================== +// Sum1D implementation +//============================================================================== + +Sum1D::Sum1D(hid_t group) +{ + // Get number of functions + int n; + read_attribute(group, "n", n); + + // Get each function + for (int i = 0; i < n; ++i) { + auto dset_name = fmt::format("func_{}", i + 1); + functions_.push_back(read_function(group, dset_name.c_str())); + } +} + +double Sum1D::operator()(double x) const +{ + double result = 0.0; + for (auto& func : functions_) { + result += (*func)(x); + } + return result; +} + } // namespace openmc diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index d53ad3879..27b750b93 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -118,6 +118,12 @@ void close_group(hid_t group_id) fatal_error("Failed to close group"); } +void close_object(hid_t obj_id) +{ + if (H5Oclose(obj_id) < 0) + fatal_error("Failed to close object"); +} + int dataset_ndims(hid_t dset) { hid_t dspace = H5Dget_space(dset); @@ -394,6 +400,12 @@ hid_t open_group(hid_t group_id, const char* name) return H5Gopen(group_id, name, H5P_DEFAULT); } +hid_t open_object(hid_t group_id, const std::string& name) +{ + ensure_exists(group_id, name.c_str()); + return H5Oopen(group_id, name.c_str(), H5P_DEFAULT); +} + void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer) { hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT); diff --git a/src/mesh.cpp b/src/mesh.cpp index e8f79aa37..55fbc8e29 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -220,21 +220,58 @@ void UnstructuredMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "type", mesh_type); write_dataset(mesh_group, "filename", filename_); write_dataset(mesh_group, "library", this->library()); - // write volume of each element - vector tet_vols; - xt::xtensor centroids({static_cast(this->n_bins()), 3}); - for (int i = 0; i < this->n_bins(); i++) { - tet_vols.emplace_back(this->volume(i)); - auto c = this->centroid(i); - xt::view(centroids, i, xt::all()) = xt::xarray({c.x, c.y, c.z}); - } - - write_dataset(mesh_group, "volumes", tet_vols); - write_dataset(mesh_group, "centroids", centroids); if (specified_length_multiplier_) write_dataset(mesh_group, "length_multiplier", length_multiplier_); + // write vertex coordinates + xt::xtensor vertices({static_cast(this->n_vertices()), 3}); + for (int i = 0; i < this->n_vertices(); i++) { + auto v = this->vertex(i); + xt::view(vertices, i, xt::all()) = xt::xarray({v.x, v.y, v.z}); + } + write_dataset(mesh_group, "vertices", vertices); + + int num_elem_skipped = 0; + + // write element types and connectivity + vector volumes; + xt::xtensor connectivity ({static_cast(this->n_bins()), 8}); + xt::xtensor elem_types ({static_cast(this->n_bins()), 1}); + for (int i = 0; i < this->n_bins(); i++) { + auto conn = this->connectivity(i); + + volumes.emplace_back(this->volume(i)); + + // write linear tet element + if (conn.size() == 4) { + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_TET); + xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], + -1, -1, -1, -1}); + // write linear hex element + } else if (conn.size() == 8) { + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_HEX); + xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], + conn[4], conn[5], conn[6], conn[7]}); + } else { + num_elem_skipped++; + xt::view(elem_types, i, xt::all()) = static_cast(ElementType::UNSUPPORTED); + xt::view(connectivity, i, xt::all()) = -1; + } + } + + // warn users that some elements were skipped + if (num_elem_skipped > 0) { + warning(fmt::format("The connectivity of {} elements " + "on mesh {} were not written " + "because they are not of type linear tet/hex.", + num_elem_skipped, this->id_)); + } + + write_dataset(mesh_group, "volumes", volumes); + write_dataset(mesh_group, "connectivity", connectivity); + write_dataset(mesh_group, "element_types", elem_types); + close_group(mesh_group); } @@ -1775,6 +1812,14 @@ void MOABMesh::initialize() filename_); } + // set member range of vertices + int vertex_dim = 0; + rval = mbi_->get_entities_by_dimension(0, vertex_dim, verts_); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get all vertex handles"); + } + + // make an entity set for all tetrahedra // this is used for convenience later in output rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_); @@ -2124,6 +2169,14 @@ std::pair, vector> MOABMesh::plot( return {}; } +int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const { + int idx = vert - verts_[0]; + if (idx >= n_vertices()) { + fatal_error(fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices())); + } + return idx; +} + int MOABMesh::get_bin_from_ent_handle(moab::EntityHandle eh) const { int bin = eh - ehs_[0]; @@ -2191,6 +2244,46 @@ Position MOABMesh::centroid(int bin) const return {centroid[0], centroid[1], centroid[2]}; } +int MOABMesh::n_vertices() const { + return verts_.size(); +} + +Position MOABMesh::vertex(int id) const { + + moab::ErrorCode rval; + + moab::EntityHandle vert = verts_[id]; + + moab::CartVect coords; + rval = mbi_->get_coords(&vert, 1, coords.array()); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get the coordinates of a vertex."); + } + + return {coords[0], coords[1], coords[2]}; +} + +std::vector MOABMesh::connectivity(int bin) const { + moab::ErrorCode rval; + + auto tet = get_ent_handle_from_bin(bin); + + // look up the tet connectivity + vector conn; + rval = mbi_->get_connectivity(&tet, 1, conn); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get connectivity of a mesh element."); + return {}; + } + + std::vector verts(4); + for (int i = 0; i < verts.size(); i++) { + verts[i] = get_vert_idx_from_handle(conn[i]); + } + + return verts; +} + std::pair MOABMesh::get_score_tags( std::string score) const { @@ -2324,16 +2417,37 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { + // filename_ and length_multiplier_ will already be set by the UnstructuredMesh constructor + set_mesh_pointer_from_filename(filename_); + set_length_multiplier(length_multiplier_); initialize(); } -LibMesh::LibMesh(const std::string& filename, double length_multiplier) +// create the mesh from a pointer to a libMesh Mesh +LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier) { - filename_ = filename; + m_ = &input_mesh; set_length_multiplier(length_multiplier); initialize(); } +// create the mesh from an input file +LibMesh::LibMesh(const std::string& filename, double length_multiplier) +{ + set_mesh_pointer_from_filename(filename); + set_length_multiplier(length_multiplier); + initialize(); +} + +void LibMesh::set_mesh_pointer_from_filename(const std::string& filename) +{ + filename_ = filename; + unique_m_ = make_unique(*settings::libmesh_comm, n_dimension_); + m_ = unique_m_.get(); + m_->read(filename_); +} + +// intialize from mesh file void LibMesh::initialize() { if (!settings::libmesh_comm) { @@ -2344,14 +2458,14 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; - m_ = make_unique(*settings::libmesh_comm, n_dimension_); - m_->read(filename_); - if (specified_length_multiplier_) { libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); } - - m_->prepare_for_use(); + // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh. + // Otherwise assume that it is prepared by its owning application + if (unique_m_) { + m_->prepare_for_use(); + } // ensure that the loaded mesh is 3 dimensional if (m_->mesh_dimension() != n_dimension_) { @@ -2394,6 +2508,27 @@ Position LibMesh::centroid(int bin) const return {centroid(0), centroid(1), centroid(2)}; } +int LibMesh::n_vertices() const +{ + return m_->n_nodes(); +} + +Position LibMesh::vertex(int vertex_id) const +{ + const auto node_ref = m_->node_ref(vertex_id); + return {node_ref(0), node_ref(1), node_ref(2)}; +} + +std::vector LibMesh::connectivity(int elem_id) const +{ + std::vector conn; + const auto* elem_ptr = m_->elem_ptr(elem_id); + for (int i = 0; i < elem_ptr->n_nodes(); i++) { + conn.push_back(elem_ptr->node_id(i)); + } + return conn; +} + std::string LibMesh::library() const { return mesh_lib_type; diff --git a/src/output.cpp b/src/output.cpp index ece475de0..56fea4db9 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -2,6 +2,7 @@ #include // for transform, max #include // for strlen +#include // for stdout #include // for time, localtime #include #include // for setw, setprecision, put_time @@ -93,7 +94,8 @@ void title() // Write number of OpenMP threads fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); #endif - std::cout << std::endl; + fmt::print("\n"); + std::fflush(stdout); } //============================================================================== @@ -133,7 +135,8 @@ void header(const char* msg, int level) // Print header based on verbosity level. if (settings::verbosity >= level) - std::cout << '\n' << out << "\n" << std::endl; + fmt::print("\n{}\n\n", out); + std::fflush(stdout); } //============================================================================== @@ -379,7 +382,8 @@ void print_generation() if (n > 1) { fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std); } - std::cout << std::endl; + fmt::print("\n"); + std::fflush(stdout); } //============================================================================== @@ -546,7 +550,8 @@ void print_results() fmt::print(" Leakage Fraction = {:.5f}\n", gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n); } - std::cout << std::endl; + fmt::print("\n"); + std::fflush(stdout); } //============================================================================== diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 3677f6ed3..0b8e1ab42 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -332,4 +332,40 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } +//============================================================================== +// MixedElasticAE implementation +//============================================================================== + +MixedElasticAE::MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs) + : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs) +{ + // Read incoherent elastic distribution + hid_t incoherent_group = open_group(group, "incoherent"); + std::string temp; + read_attribute(incoherent_group, "type", temp); + if (temp == "incoherent_elastic") { + incoherent_dist_ = make_unique(incoherent_group); + } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(&incoh_xs); + incoherent_dist_ = + make_unique(incoherent_group, xs->x()); + } + close_group(incoherent_group); +} + +void MixedElasticAE::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // Evaluate coherent and incoherent elastic cross sections + double xs_coh = coherent_xs_(E_in); + double xs_incoh = incoherent_xs_(E_in); + + if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) { + coherent_dist_.sample(E_in, E_out, mu, seed); + } else { + incoherent_dist_->sample(E_in, E_out, mu, seed); + } +} + } // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index 11ede5fc4..3716720df 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -399,4 +399,23 @@ void free_memory_source() model::external_sources.clear(); } +//============================================================================== +// C API +//============================================================================== + +extern "C" int openmc_sample_external_source( + size_t n, uint64_t* seed, void* sites) +{ + 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; +} + } // namespace openmc diff --git a/src/state_point.cpp b/src/state_point.cpp index 9ea045915..648d1a249 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -818,6 +818,16 @@ void write_unstructured_mesh_results() if (!umesh->output_) continue; + if (umesh->library() == "moab") { + if (mpi::master) + warning(fmt::format( + "Output for a MOAB mesh (mesh {}) was " + "requested but will not be written. Please use the Python " + "API to generated the desired VTK tetrahedral mesh.", + umesh->id_)); + continue; + } + // if this tally has more than one filter, print // warning and skip writing the mesh if (tally->filters().size() > 1) { @@ -889,12 +899,10 @@ void write_unstructured_mesh_results() std::string filename = fmt::format("tally_{0}.{1:0{2}}", tally->id_, simulation::current_batch, batch_width); - if (umesh->library() == "moab" && !mpi::master) - continue; - // Write the unstructured mesh and data to file umesh->write(filename); + // remove score data added for this mesh write umesh->remove_scores(); } } diff --git a/src/thermal.cpp b/src/thermal.cpp index 1101d5485..1a9592d0a 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -210,14 +210,22 @@ ThermalData::ThermalData(hid_t group) if (temp == "coherent_elastic") { auto xs = dynamic_cast(elastic_.xs.get()); elastic_.distribution = make_unique(*xs); - } else { - if (temp == "incoherent_elastic") { - elastic_.distribution = make_unique(dgroup); - } else if (temp == "incoherent_elastic_discrete") { - auto xs = dynamic_cast(elastic_.xs.get()); - elastic_.distribution = - make_unique(dgroup, xs->x()); - } + } else if (temp == "incoherent_elastic") { + elastic_.distribution = make_unique(dgroup); + } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(elastic_.xs.get()); + elastic_.distribution = + make_unique(dgroup, xs->x()); + } else if (temp == "mixed_elastic") { + // Get coherent/incoherent cross sections + auto mixed_xs = dynamic_cast(elastic_.xs.get()); + const auto& coh_xs = + dynamic_cast(mixed_xs->functions(0).get()); + const auto& incoh_xs = mixed_xs->functions(1).get(); + + // Create mixed elastic distribution + elastic_.distribution = + make_unique(dgroup, *coh_xs, *incoh_xs); } close_group(elastic_group); diff --git a/tests/regression_tests/dagmc/universes/inputs_true.dat b/tests/regression_tests/dagmc/universes/inputs_true.dat index 4443f9a35..2165a78b6 100644 --- a/tests/regression_tests/dagmc/universes/inputs_true.dat +++ b/tests/regression_tests/dagmc/universes/inputs_true.dat @@ -1,8 +1,8 @@ - + - + 24.0 24.0 2 2 -24.0 -24.0 @@ -10,12 +10,12 @@ 9 9 9 9 - - - - - - + + + + + + diff --git a/tests/regression_tests/dagmc/universes/test.py b/tests/regression_tests/dagmc/universes/test.py index 6726b4b44..f2eb3882f 100644 --- a/tests/regression_tests/dagmc/universes/test.py +++ b/tests/regression_tests/dagmc/universes/test.py @@ -46,11 +46,32 @@ class DAGMCUniverseTest(PyAPITestHarness): # create the DAGMC universe pincell_univ = openmc.DAGMCUniverse(filename='dagmc.h5m', auto_geom_ids=True) + # creates another DAGMC universe, this time with within a bounded cell + bound_pincell_universe = openmc.DAGMCUniverse(filename='dagmc.h5m').bounded_universe() + # uses the bound_dag_cell as the root argument to test the type checks in openmc.Geometry + bound_pincell_geometry = openmc.Geometry(root=bound_pincell_universe) + # assigns the bound_dag_geometry to the model to test the type checks in model.Geometry setter + model.Geometry = bound_pincell_geometry + # checks that the bounding box is calculated correctly bounding_box = pincell_univ.bounding_box assert bounding_box[0].tolist() == [-25., -25., -25.] assert bounding_box[1].tolist() == [25., 25., 25.] + # checks that the bounding region is six surfaces each with a vacuum boundary type + b_region = pincell_univ.bounding_region(bounded_type='box', boundary_type='vacuum') + assert isinstance(b_region, openmc.Region) + assert len(b_region.get_surfaces()) == 6 + for surface in list(b_region.get_surfaces().values()): + assert surface.boundary_type == 'vacuum' + + # checks that the bounding region is a single surface with a reflective boundary type + b_region = pincell_univ.bounding_region(bounded_type='sphere', boundary_type='reflective') + assert isinstance(b_region, openmc.Region) + assert len(b_region.get_surfaces()) == 1 + for surface in list(b_region.get_surfaces().values()): + assert surface.boundary_type == 'reflective' + # create a 2 x 2 lattice using the DAGMC pincell pitch = np.asarray((24.0, 24.0)) lattice = openmc.RectLattice() diff --git a/tests/regression_tests/unstructured_mesh/inputs_true.dat b/tests/regression_tests/unstructured_mesh/inputs_true.dat new file mode 100644 index 000000000..e556782ad --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/inputs_true.dat @@ -0,0 +1,91 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_hexes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true0.dat b/tests/regression_tests/unstructured_mesh/inputs_true0.dat index 1ed60d507..2e2594cde 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true0.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true0.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true1.dat b/tests/regression_tests/unstructured_mesh/inputs_true1.dat index 67dde56a9..84a3b186b 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true1.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true1.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true10.dat b/tests/regression_tests/unstructured_mesh/inputs_true10.dat index 9ca47a356..e9272a333 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true10.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true10.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true11.dat b/tests/regression_tests/unstructured_mesh/inputs_true11.dat index 683e1fade..0b89d280c 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true11.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true11.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true12.dat b/tests/regression_tests/unstructured_mesh/inputs_true12.dat index 8ce9f3cf2..b3673a254 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true12.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true12.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true13.dat b/tests/regression_tests/unstructured_mesh/inputs_true13.dat index bca0bee4c..c466396b6 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true13.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true13.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true14.dat b/tests/regression_tests/unstructured_mesh/inputs_true14.dat index 226331ba8..8c6a5120a 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true14.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true14.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true15.dat b/tests/regression_tests/unstructured_mesh/inputs_true15.dat index a6a084165..1b9abbd83 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true15.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true15.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true2.dat b/tests/regression_tests/unstructured_mesh/inputs_true2.dat index 4b442c757..7136d485a 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true2.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true2.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true3.dat b/tests/regression_tests/unstructured_mesh/inputs_true3.dat index 52e53498e..23900a060 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true3.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true3.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true8.dat b/tests/regression_tests/unstructured_mesh/inputs_true8.dat index e484c95a2..a0cbbcae7 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true8.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true8.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/inputs_true9.dat b/tests/regression_tests/unstructured_mesh/inputs_true9.dat index 5e83d71bd..930bd8579 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true9.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true9.dat @@ -15,12 +15,12 @@ - - - - - - + + + + + + diff --git a/tests/regression_tests/unstructured_mesh/test.py b/tests/regression_tests/unstructured_mesh/test.py index daebe7948..090a47599 100644 --- a/tests/regression_tests/unstructured_mesh/test.py +++ b/tests/regression_tests/unstructured_mesh/test.py @@ -1,6 +1,8 @@ +import filecmp import glob from itertools import product import os +import warnings import openmc import openmc.lib @@ -9,15 +11,29 @@ import numpy as np import pytest from tests.testing_harness import PyAPITestHarness -TETS_PER_VOXEL = 12 - class UnstructuredMeshTest(PyAPITestHarness): - def __init__(self, statepoint_name, model, inputs_true, holes): + ELEM_PER_VOXEL = 12 + + def __init__(self, + statepoint_name, + model, + inputs_true='inputs_true.dat', + holes=False, + scale_factor=10.0): super().__init__(statepoint_name, model, inputs_true) self.holes = holes # holes in the test mesh + self.scale_bounding_cell(scale_factor) + + def scale_bounding_cell(self, scale_factor): + geometry = self._model.geometry + for surface in geometry.get_all_surfaces().values(): + if surface.boundary_type != 'vacuum': + continue + for coeff in surface._coefficients: + surface._coefficients[coeff] *= scale_factor def _compare_results(self): with openmc.StatePoint(self._sp_name) as sp: @@ -28,32 +44,27 @@ class UnstructuredMeshTest(PyAPITestHarness): flt = tally.find_filter(openmc.MeshFilter) if isinstance(flt.mesh, openmc.RegularMesh): - reg_mesh_data, reg_mesh_std_dev = self.get_mesh_tally_data(tally) + reg_mesh_data = self.get_mesh_tally_data(tally) if self.holes: reg_mesh_data = np.delete(reg_mesh_data, self.holes) - reg_mesh_std_dev = np.delete(reg_mesh_std_dev, self.holes) else: umesh_tally = tally - unstructured_data, unstructured_std_dev = self.get_mesh_tally_data(tally, True) + unstructured_data = self.get_mesh_tally_data(tally, True) - # we expect these results to be the same to within at least ten - # decimal places - decimals = 10 if umesh_tally.estimator == 'collision' else 8 - np.testing.assert_array_almost_equal(unstructured_data, - reg_mesh_data, - decimals) + # we expect these results to be the same to within at least ten + # decimal places + decimals = 10 if umesh_tally.estimator == 'collision' else 8 + np.testing.assert_array_almost_equal(np.sort(unstructured_data), + np.sort(reg_mesh_data), + decimals) - @staticmethod - def get_mesh_tally_data(tally, structured=False): + def get_mesh_tally_data(self, tally, structured=False): data = tally.get_reshaped_data(value='mean') - std_dev = tally.get_reshaped_data(value='std_dev') if structured: - data.shape = (data.size // TETS_PER_VOXEL, TETS_PER_VOXEL) - std_dev.shape = (std_dev.size // TETS_PER_VOXEL, TETS_PER_VOXEL) + data = data.reshape((-1, self.ELEM_PER_VOXEL)) else: data.shape = (data.size, 1) - std_dev.shape = (std_dev.size, 1) - return np.sum(data, axis=1), np.sum(std_dev, axis=1) + return np.sum(data, axis=1) def _cleanup(self): super()._cleanup() @@ -64,35 +75,11 @@ class UnstructuredMeshTest(PyAPITestHarness): os.remove(f) -param_values = (['libmesh', 'moab'], # mesh libraries - ['collision', 'tracklength'], # estimators - [True, False], # geometry outside of the mesh - [(333, 90, 77), None]) # location of holes in the mesh -test_cases = [] -for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)): - test_cases.append({'library' : lib, - 'estimator' : estimator, - 'external_geom' : ext_geom, - 'holes' : holes, - 'inputs_true' : 'inputs_true{}.dat'.format(i)}) - - -@pytest.mark.parametrize("test_opts", test_cases) -def test_unstructured_mesh(test_opts): - +@pytest.fixture +def model(): openmc.reset_auto_ids() - # skip the test if the library is not enabled - if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled(): - pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.") - - if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): - pytest.skip("LibMesh is not enabled in this build.") - - # skip the tracklength test for libmesh - if test_opts['library'] == 'libmesh' and \ - test_opts['estimator'] == 'tracklength': - pytest.skip("Tracklength tallies are not supported using libmesh.") + model = openmc.Model() ### Materials ### materials = openmc.Materials() @@ -113,7 +100,7 @@ def test_unstructured_mesh(test_opts): water_mat.set_density("atom/b-cm", 0.07416) materials.append(water_mat) - materials.export_to_xml() + model.materials = materials ### Geometry ### fuel_min_x = openmc.XPlane(-5.0, name="minimum x") @@ -149,29 +136,26 @@ def test_unstructured_mesh(test_opts): +clad_min_z & -clad_max_z) clad_cell.fill = zirc_mat - if test_opts['external_geom']: - bounds = (15, 15, 15) - else: - bounds = (10, 10, 10) - - water_min_x = openmc.XPlane(x0=-bounds[0], + # set bounding cell dimension to one + # this will be updated later according to the test case parameters + water_min_x = openmc.XPlane(x0=-1.0, name="minimum x", boundary_type='vacuum') - water_max_x = openmc.XPlane(x0=bounds[0], + water_max_x = openmc.XPlane(x0=1.0, name="maximum x", boundary_type='vacuum') - water_min_y = openmc.YPlane(y0=-bounds[1], + water_min_y = openmc.YPlane(y0=-1.0, name="minimum y", boundary_type='vacuum') - water_max_y = openmc.YPlane(y0=bounds[1], + water_max_y = openmc.YPlane(y0=1.0, name="maximum y", boundary_type='vacuum') - water_min_z = openmc.ZPlane(z0=-bounds[2], + water_min_z = openmc.ZPlane(z0=-1.0, name="minimum z", boundary_type='vacuum') - water_max_z = openmc.ZPlane(z0=bounds[2], + water_max_z = openmc.ZPlane(z0=1.0, name="maximum z", boundary_type='vacuum') @@ -185,9 +169,9 @@ def test_unstructured_mesh(test_opts): water_cell.fill = water_mat # create a containing universe - geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell]) + model.geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell]) - ### Tallies ### + ### Reference Tally ### # create meshes and mesh filters regular_mesh = openmc.RegularMesh() @@ -196,29 +180,11 @@ def test_unstructured_mesh(test_opts): regular_mesh.upper_right = (10.0, 10.0, 10.0) regular_mesh_filter = openmc.MeshFilter(mesh=regular_mesh) - - if test_opts['holes']: - mesh_filename = "test_mesh_tets_w_holes.e" - else: - mesh_filename = "test_mesh_tets.e" - - uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library']) - uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) - - # create tallies - tallies = openmc.Tallies() - regular_mesh_tally = openmc.Tally(name="regular mesh tally") regular_mesh_tally.filters = [regular_mesh_filter] regular_mesh_tally.scores = ['flux'] - regular_mesh_tally.estimator = test_opts['estimator'] - tallies.append(regular_mesh_tally) - uscd_tally = openmc.Tally(name="unstructured mesh tally") - uscd_tally.filters = [uscd_filter] - uscd_tally.scores = ['flux'] - uscd_tally.estimator = test_opts['estimator'] - tallies.append(uscd_tally) + model.tallies = openmc.Tallies([regular_mesh_tally]) ### Settings ### settings = openmc.Settings() @@ -236,13 +202,92 @@ def test_unstructured_mesh(test_opts): source = openmc.Source(space=space, energy=energy) settings.source = source - model = openmc.model.Model(geometry=geometry, - materials=materials, - tallies=tallies, - settings=settings) + model.settings = settings + + return model + + +param_values = (['libmesh', 'moab'], # mesh libraries + ['collision', 'tracklength'], # estimators + [True, False], # geometry outside of the mesh + [(333, 90, 77), None]) # location of holes in the mesh +test_cases = [] +for i, (lib, estimator, ext_geom, holes) in enumerate(product(*param_values)): + test_cases.append({'library' : lib, + 'estimator' : estimator, + 'external_geom' : ext_geom, + 'holes' : holes, + 'inputs_true' : 'inputs_true{}.dat'.format(i)}) + + +@pytest.mark.parametrize("test_opts", test_cases) +def test_unstructured_mesh_tets(model, test_opts): + # skip the test if the library is not enabled + if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip("DAGMC (and MOAB) mesh not enbaled in this build.") + + if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip("LibMesh is not enabled in this build.") + + # skip the tracklength test for libmesh + if test_opts['library'] == 'libmesh' and \ + test_opts['estimator'] == 'tracklength': + pytest.skip("Tracklength tallies are not supported using libmesh.") + + if test_opts['holes']: + mesh_filename = "test_mesh_tets_w_holes.e" + else: + mesh_filename = "test_mesh_tets.e" + + # add reference mesh tally + regular_mesh_tally = model.tallies[0] + regular_mesh_tally.estimator = test_opts['estimator'] + + # add analagous unstructured mesh tally + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_opts['library']) + uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) + + # create tallies + uscd_tally = openmc.Tally(name="unstructured mesh tally") + uscd_tally.filters = [uscd_filter] + uscd_tally.scores = ['flux'] + uscd_tally.estimator = test_opts['estimator'] + model.tallies.append(uscd_tally) + + # modify model geometry according to test opts + if test_opts['external_geom']: + scale_factor = 15.0 + else: + scale_factor = 10.0 harness = UnstructuredMeshTest('statepoint.10.h5', model, test_opts['inputs_true'], - test_opts['holes']) + test_opts['holes'], + scale_factor) harness.main() + + +@pytest.mark.skipif(not openmc.lib._libmesh_enabled(), + reason='LibMesh is not enabled in this build.') +def test_unstructured_mesh_hexes(model): + regular_mesh_tally = model.tallies[0] + regular_mesh_tally.estimator = 'collision' + + # add analagous unstructured mesh tally + uscd_mesh = openmc.UnstructuredMesh('test_mesh_hexes.e', 'libmesh') + uscd_filter = openmc.MeshFilter(mesh=uscd_mesh) + + # create tallies + uscd_tally = openmc.Tally(name="unstructured mesh tally") + uscd_tally.filters = [uscd_filter] + uscd_tally.scores = ['flux'] + uscd_tally.estimator = 'collision' + model.tallies.append(uscd_tally) + + harness = UnstructuredMeshTest('statepoint.10.h5', + model) + harness.ELEM_PER_VOXEL = 1 + + harness.main() + diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e new file mode 120000 index 000000000..421bf6a89 --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.e @@ -0,0 +1 @@ +test_mesh_hexes.exo \ No newline at end of file diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo new file mode 100644 index 000000000..03682c275 Binary files /dev/null and b/tests/regression_tests/unstructured_mesh/test_mesh_hexes.exo differ diff --git a/tests/unit_tests/mesh_to_vtk/__init__.py b/tests/unit_tests/mesh_to_vtk/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/unit_tests/mesh_to_vtk/hexes.exo b/tests/unit_tests/mesh_to_vtk/hexes.exo new file mode 100644 index 000000000..c16138e90 Binary files /dev/null and b/tests/unit_tests/mesh_to_vtk/hexes.exo differ diff --git a/tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk b/tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk new file mode 100644 index 000000000..041d7065d --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/libmesh_hexes_ref.vtk @@ -0,0 +1,76 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 54 double +-1 -1 2.5 -1 -1 1.5 -1 0 1.5 +-1 0 2.5 0 -1 2.5 0 -1 1.5 +0 0 1.5 0 0 2.5 -1 -1 0.5 +-1 0 0.5 0 -1 0.5 0 0 0.5 +-1 -1 -0.5 -1 0 -0.5 0 -1 -0.5 +0 0 -0.5 -1 -1 -1.5 -1 0 -1.5 +0 -1 -1.5 0 0 -1.5 -1 -1 -2.5 +-1 0 -2.5 0 -1 -2.5 0 0 -2.5 +-1 1 1.5 -1 1 2.5 0 1 1.5 +0 1 2.5 -1 1 0.5 0 1 0.5 +-1 1 -0.5 0 1 -0.5 -1 1 -1.5 +0 1 -1.5 -1 1 -2.5 0 1 -2.5 +1 -1 2.5 1 -1 1.5 1 0 1.5 +1 0 2.5 1 -1 0.5 1 0 0.5 +1 -1 -0.5 1 0 -0.5 1 -1 -1.5 +1 0 -1.5 1 -1 -2.5 1 0 -2.5 +1 1 1.5 1 1 2.5 1 1 0.5 +1 1 -0.5 1 1 -1.5 1 1 -2.5 + +CELLS 21 160 +OFFSETS vtktypeint64 +0 8 16 24 32 40 48 56 64 +72 80 88 96 104 112 120 128 136 +144 152 160 +CONNECTIVITY vtktypeint64 +0 1 2 3 4 5 6 7 1 +8 9 2 5 10 11 6 8 12 +13 9 10 14 15 11 12 16 17 +13 14 18 19 15 16 20 21 17 +18 22 23 19 3 2 24 25 7 +6 26 27 2 9 28 24 6 11 +29 26 9 13 30 28 11 15 31 +29 13 17 32 30 15 19 33 31 +17 21 34 32 19 23 35 33 4 +5 6 7 36 37 38 39 5 10 +11 6 37 40 41 38 10 14 15 +11 40 42 43 41 14 18 19 15 +42 44 45 43 18 22 23 19 44 +46 47 45 7 6 26 27 39 38 +48 49 6 11 29 26 38 41 50 +48 11 15 31 29 41 43 51 50 +15 19 33 31 43 45 52 51 19 +23 35 33 45 47 53 52 +CELL_TYPES 20 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 +12 + +CELL_DATA 20 +FIELD FieldData 1 +ids 1 20 double +0 1 2 3 4 5 6 7 8 +9 10 11 12 13 14 15 16 17 +18 19 diff --git a/tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk b/tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk new file mode 100644 index 000000000..90f226121 --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/libmesh_tets_ref.vtk @@ -0,0 +1,292 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 58 double +-0.02593964576 -1 -1.1195739536 -0.40239958217 -0.40962166746 -1.9256035383 -1 0.02593964576 -1.1195739536 +-0.2248833639 0.23189144433 -1.2841109811 0.02593964576 1 -1.1195739536 1 -0.02593964576 -1.1195739536 +-0.042928712675 0.066425810853 -0.54270728017 -1 1 -1.5 -1 1 -0.5 +-0.49265312381 0.49332495053 -1.9656359509 -0.02608137991 -1 1.121521147 0.30281888252 -0.30293622513 1.8067963894 +1 -0.02608137991 1.121521147 -0.23236342389 0.26314805583 1.2867556783 0.074337770194 0.074337770194 2.5 +-0.40005770996 -0.36909660956 1.8768871762 -0.4896743493 0.49451064621 1.952829049 -1 1 1.5 +0.02608137991 1 1.121521147 -1 1 0.5 -0.038567136385 0.14628887118 0.55111724874 +-1 0.02608137991 1.121521147 0 -1 -2.5 0.37543954218 -0.36967061971 -1.8691105401 +1 -1 -2.5 0.074337770194 -0.074337770194 -2.5 1 0 -2.5 +0.4135797166 0.38351746262 -1.9274856063 1 -1 -1.5 -1 0 -2.5 +-1 -1 -1.5 -1 -1 -2.5 -1 0 2.5 +-1 1 2.5 0 1 2.5 0.46519182209 0.49356920409 1.95286052 +0 -1 2.5 -1 -1 1.5 -1 -1 2.5 +1 -1 2.5 1 -1 1.5 1 0 2.5 +1 1 2.5 1 1 1.5 0 1 -2.5 +-1 1 -2.5 1 -1 -0.5 1 1 -0.5 +0.038669824604 1 0.0046021677516 -1 -1 0.5 -0.038669824604 -1 0.0046021677516 +1 1 0.5 -1 0.038669824604 0.0046021677516 1 -0.038669824604 0.0046021677516 +1 1 -2.5 1 1 -1.5 1 -1 0.5 +-1 -1 -0.5 +CELLS 155 616 +OFFSETS vtktypeint64 +0 4 8 12 16 20 24 28 32 +36 40 44 48 52 56 60 64 68 +72 76 80 84 88 92 96 100 104 +108 112 116 120 124 128 132 136 140 +144 148 152 156 160 164 168 172 176 +180 184 188 192 196 200 204 208 212 +216 220 224 228 232 236 240 244 248 +252 256 260 264 268 272 276 280 284 +288 292 296 300 304 308 312 316 320 +324 328 332 336 340 344 348 352 356 +360 364 368 372 376 380 384 388 392 +396 400 404 408 412 416 420 424 428 +432 436 440 444 448 452 456 460 464 +468 472 476 480 484 488 492 496 500 +504 508 512 516 520 524 528 532 536 +540 544 548 552 556 560 564 568 572 +576 580 584 588 592 596 600 604 608 +612 616 +CONNECTIVITY vtktypeint64 +0 1 2 3 4 3 5 6 7 +4 8 3 7 3 2 9 8 3 +4 6 10 11 12 13 14 13 15 +16 17 13 18 16 19 18 13 20 +10 13 21 15 22 1 0 23 22 +24 25 23 25 1 3 9 26 23 +5 27 24 28 26 23 0 1 3 +23 29 30 31 1 22 31 30 1 +2 3 1 9 29 2 1 9 22 +0 28 23 32 15 21 16 14 13 +11 15 33 17 34 16 33 34 32 +16 34 16 18 35 10 11 13 15 +18 13 12 35 36 11 10 15 36 +37 38 15 32 38 37 15 32 21 +17 16 36 10 37 15 36 39 40 +11 41 40 39 11 41 12 11 35 +17 21 13 16 10 13 12 20 12 +13 11 35 14 13 16 35 34 17 +18 16 42 34 43 35 41 12 40 +11 14 41 11 35 42 41 14 35 +0 3 5 23 0 5 3 6 44 +4 9 27 29 2 30 1 44 45 +29 9 29 45 7 9 0 46 5 +6 47 48 4 6 0 3 2 6 +4 5 3 27 49 10 50 20 8 +2 3 6 51 18 48 20 49 50 +52 20 19 13 21 20 51 53 12 +20 19 48 18 20 25 3 1 23 +7 4 3 9 25 3 23 27 44 +4 7 9 25 23 26 27 54 55 +44 27 54 26 55 27 18 12 13 +20 14 16 34 35 19 21 52 20 +30 2 0 1 22 25 29 1 5 +23 3 27 25 26 44 27 22 29 +31 1 54 44 26 27 25 44 9 +27 22 28 24 23 26 5 55 27 +25 24 26 23 22 25 1 23 25 +44 29 9 25 29 1 9 48 6 +53 20 8 4 48 6 10 40 56 +12 10 40 12 11 17 19 18 13 +30 57 0 2 25 9 3 27 7 +8 2 3 26 28 5 23 29 7 +2 9 28 0 5 23 37 10 49 +21 49 50 57 52 56 12 53 20 +10 12 56 20 43 18 12 35 32 +37 21 15 14 11 13 35 19 52 +8 48 19 52 48 20 41 39 36 +11 50 6 52 20 46 53 5 6 +4 3 9 27 44 7 45 9 55 +47 4 5 28 0 46 5 8 52 +2 6 57 2 52 6 50 56 53 +20 10 56 50 20 57 50 0 6 +37 10 21 15 33 32 17 16 34 +18 43 35 57 52 50 6 51 48 +47 53 18 16 13 35 50 56 46 +53 14 36 32 15 14 11 36 15 +17 21 19 13 10 21 13 20 51 +12 18 20 50 53 46 6 14 15 +32 16 42 43 41 35 14 41 36 +11 49 52 21 20 36 40 10 11 +47 53 48 6 49 21 10 20 21 +15 13 16 51 48 53 20 55 5 +4 27 44 55 4 27 47 4 5 +6 8 48 52 6 47 5 53 6 +50 53 6 20 57 0 2 6 50 +46 0 6 52 6 48 20 42 14 +34 35 43 18 51 12 41 43 12 +35 22 30 0 1 14 32 34 16 +36 38 32 15 +CELL_TYPES 154 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 + +CELL_DATA 154 +FIELD FieldData 1 +ids 1 154 double +0 1 2 3 4 5 6 7 8 +9 10 11 12 13 14 15 16 17 +18 19 20 21 22 23 24 25 26 +27 28 29 30 31 32 33 34 35 +36 37 38 39 40 41 42 43 44 +45 46 47 48 49 50 51 52 53 +54 55 56 57 58 59 60 61 62 +63 64 65 66 67 68 69 70 71 +72 73 74 75 76 77 78 79 80 +81 82 83 84 85 86 87 88 89 +90 91 92 93 94 95 96 97 98 +99 100 101 102 103 104 105 106 107 +108 109 110 111 112 113 114 115 116 +117 118 119 120 121 122 123 124 125 +126 127 128 129 130 131 132 133 134 +135 136 137 138 139 140 141 142 143 +144 145 146 147 148 149 150 151 152 +153 diff --git a/tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk b/tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk new file mode 120000 index 000000000..91e356480 --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/moab_tets_ref.vtk @@ -0,0 +1 @@ +libmesh_tets_ref.vtk \ No newline at end of file diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py new file mode 100644 index 000000000..5f6cd6a41 --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -0,0 +1,69 @@ +import filecmp +from itertools import product +from pathlib import Path + +import numpy as np + +import openmc +import openmc.lib +import pytest + +pytest.importorskip('vtk') + + +def ids_func(param): + return f"{param['library']}_{param['elem_type']}" + +test_params = (['libmesh', 'moab'], + ['tets', 'hexes']) + +test_cases = [ + {'library' : library, 'elem_type' : elem_type} + for library, elem_type in product(*test_params) +] + +@pytest.mark.parametrize("test_opts", test_cases, ids=ids_func) +def test_unstructured_mesh_to_vtk(run_in_tmpdir, request, test_opts): + + if test_opts['library'] == 'moab' and test_opts['elem_type'] == 'hexes': + pytest.skip('Hexes are not supported with MOAB') + + if test_opts['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip('LibMesh is not enabled in this build.') + + if test_opts['library'] == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip('DAGMC (and MOAB) mesh not enabled in this build.') + + # pull in a simple model -- just need to create the statepoint file + openmc.reset_auto_ids() + model = openmc.examples.pwr_pin_cell() + + if test_opts['elem_type'] == 'tets': + filename = Path('tets.exo') + else: + filename = Path('hexes.exo') + + # create a basic tally using the unstructured mesh + umesh = openmc.UnstructuredMesh(request.node.path.parent / filename, + test_opts['library']) + umesh.output = False + mesh_filter = openmc.MeshFilter(umesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + tally.estimator = 'collision' + model.tallies = openmc.Tallies([tally]) + sp_file = model.run() + + # check VTK output after reading mesh from statepoint file + with openmc.StatePoint(sp_file) as sp: + umesh = sp.meshes[umesh.id] + + test_data = {'ids': np.arange(umesh.n_elements)} + umesh.write_data_to_vtk('umesh.vtk', + datasets=test_data, + volume_normalization=False) + + # compare file content with reference file + ref_file = Path(f"{test_opts['library']}_{test_opts['elem_type']}_ref.vtk") + assert filecmp.cmp('umesh.vtk', request.node.path.parent / ref_file) diff --git a/tests/unit_tests/mesh_to_vtk/tets.exo b/tests/unit_tests/mesh_to_vtk/tets.exo new file mode 100644 index 000000000..add4f7cfa Binary files /dev/null and b/tests/unit_tests/mesh_to_vtk/tets.exo differ diff --git a/tests/unit_tests/mesh_to_vtk/umesh.vtk b/tests/unit_tests/mesh_to_vtk/umesh.vtk new file mode 100644 index 000000000..90f226121 --- /dev/null +++ b/tests/unit_tests/mesh_to_vtk/umesh.vtk @@ -0,0 +1,292 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 58 double +-0.02593964576 -1 -1.1195739536 -0.40239958217 -0.40962166746 -1.9256035383 -1 0.02593964576 -1.1195739536 +-0.2248833639 0.23189144433 -1.2841109811 0.02593964576 1 -1.1195739536 1 -0.02593964576 -1.1195739536 +-0.042928712675 0.066425810853 -0.54270728017 -1 1 -1.5 -1 1 -0.5 +-0.49265312381 0.49332495053 -1.9656359509 -0.02608137991 -1 1.121521147 0.30281888252 -0.30293622513 1.8067963894 +1 -0.02608137991 1.121521147 -0.23236342389 0.26314805583 1.2867556783 0.074337770194 0.074337770194 2.5 +-0.40005770996 -0.36909660956 1.8768871762 -0.4896743493 0.49451064621 1.952829049 -1 1 1.5 +0.02608137991 1 1.121521147 -1 1 0.5 -0.038567136385 0.14628887118 0.55111724874 +-1 0.02608137991 1.121521147 0 -1 -2.5 0.37543954218 -0.36967061971 -1.8691105401 +1 -1 -2.5 0.074337770194 -0.074337770194 -2.5 1 0 -2.5 +0.4135797166 0.38351746262 -1.9274856063 1 -1 -1.5 -1 0 -2.5 +-1 -1 -1.5 -1 -1 -2.5 -1 0 2.5 +-1 1 2.5 0 1 2.5 0.46519182209 0.49356920409 1.95286052 +0 -1 2.5 -1 -1 1.5 -1 -1 2.5 +1 -1 2.5 1 -1 1.5 1 0 2.5 +1 1 2.5 1 1 1.5 0 1 -2.5 +-1 1 -2.5 1 -1 -0.5 1 1 -0.5 +0.038669824604 1 0.0046021677516 -1 -1 0.5 -0.038669824604 -1 0.0046021677516 +1 1 0.5 -1 0.038669824604 0.0046021677516 1 -0.038669824604 0.0046021677516 +1 1 -2.5 1 1 -1.5 1 -1 0.5 +-1 -1 -0.5 +CELLS 155 616 +OFFSETS vtktypeint64 +0 4 8 12 16 20 24 28 32 +36 40 44 48 52 56 60 64 68 +72 76 80 84 88 92 96 100 104 +108 112 116 120 124 128 132 136 140 +144 148 152 156 160 164 168 172 176 +180 184 188 192 196 200 204 208 212 +216 220 224 228 232 236 240 244 248 +252 256 260 264 268 272 276 280 284 +288 292 296 300 304 308 312 316 320 +324 328 332 336 340 344 348 352 356 +360 364 368 372 376 380 384 388 392 +396 400 404 408 412 416 420 424 428 +432 436 440 444 448 452 456 460 464 +468 472 476 480 484 488 492 496 500 +504 508 512 516 520 524 528 532 536 +540 544 548 552 556 560 564 568 572 +576 580 584 588 592 596 600 604 608 +612 616 +CONNECTIVITY vtktypeint64 +0 1 2 3 4 3 5 6 7 +4 8 3 7 3 2 9 8 3 +4 6 10 11 12 13 14 13 15 +16 17 13 18 16 19 18 13 20 +10 13 21 15 22 1 0 23 22 +24 25 23 25 1 3 9 26 23 +5 27 24 28 26 23 0 1 3 +23 29 30 31 1 22 31 30 1 +2 3 1 9 29 2 1 9 22 +0 28 23 32 15 21 16 14 13 +11 15 33 17 34 16 33 34 32 +16 34 16 18 35 10 11 13 15 +18 13 12 35 36 11 10 15 36 +37 38 15 32 38 37 15 32 21 +17 16 36 10 37 15 36 39 40 +11 41 40 39 11 41 12 11 35 +17 21 13 16 10 13 12 20 12 +13 11 35 14 13 16 35 34 17 +18 16 42 34 43 35 41 12 40 +11 14 41 11 35 42 41 14 35 +0 3 5 23 0 5 3 6 44 +4 9 27 29 2 30 1 44 45 +29 9 29 45 7 9 0 46 5 +6 47 48 4 6 0 3 2 6 +4 5 3 27 49 10 50 20 8 +2 3 6 51 18 48 20 49 50 +52 20 19 13 21 20 51 53 12 +20 19 48 18 20 25 3 1 23 +7 4 3 9 25 3 23 27 44 +4 7 9 25 23 26 27 54 55 +44 27 54 26 55 27 18 12 13 +20 14 16 34 35 19 21 52 20 +30 2 0 1 22 25 29 1 5 +23 3 27 25 26 44 27 22 29 +31 1 54 44 26 27 25 44 9 +27 22 28 24 23 26 5 55 27 +25 24 26 23 22 25 1 23 25 +44 29 9 25 29 1 9 48 6 +53 20 8 4 48 6 10 40 56 +12 10 40 12 11 17 19 18 13 +30 57 0 2 25 9 3 27 7 +8 2 3 26 28 5 23 29 7 +2 9 28 0 5 23 37 10 49 +21 49 50 57 52 56 12 53 20 +10 12 56 20 43 18 12 35 32 +37 21 15 14 11 13 35 19 52 +8 48 19 52 48 20 41 39 36 +11 50 6 52 20 46 53 5 6 +4 3 9 27 44 7 45 9 55 +47 4 5 28 0 46 5 8 52 +2 6 57 2 52 6 50 56 53 +20 10 56 50 20 57 50 0 6 +37 10 21 15 33 32 17 16 34 +18 43 35 57 52 50 6 51 48 +47 53 18 16 13 35 50 56 46 +53 14 36 32 15 14 11 36 15 +17 21 19 13 10 21 13 20 51 +12 18 20 50 53 46 6 14 15 +32 16 42 43 41 35 14 41 36 +11 49 52 21 20 36 40 10 11 +47 53 48 6 49 21 10 20 21 +15 13 16 51 48 53 20 55 5 +4 27 44 55 4 27 47 4 5 +6 8 48 52 6 47 5 53 6 +50 53 6 20 57 0 2 6 50 +46 0 6 52 6 48 20 42 14 +34 35 43 18 51 12 41 43 12 +35 22 30 0 1 14 32 34 16 +36 38 32 15 +CELL_TYPES 154 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 + +CELL_DATA 154 +FIELD FieldData 1 +ids 1 154 double +0 1 2 3 4 5 6 7 8 +9 10 11 12 13 14 15 16 17 +18 19 20 21 22 23 24 25 26 +27 28 29 30 31 32 33 34 35 +36 37 38 39 40 41 42 43 44 +45 46 47 48 49 50 51 52 53 +54 55 56 57 58 59 60 61 62 +63 64 65 66 67 68 69 70 71 +72 73 74 75 76 77 78 79 80 +81 82 83 84 85 86 87 88 89 +90 91 92 93 94 95 96 97 98 +99 100 101 102 103 104 105 106 107 +108 109 110 111 112 113 114 115 116 +117 118 119 120 121 122 123 124 125 +126 127 128 129 130 131 132 133 134 +135 136 137 138 139 140 141 142 143 +144 145 146 147 148 149 150 151 152 +153 diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index dc4d24628..61ae61b24 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -262,3 +262,108 @@ def test_get_thermal_name(): # Names that don't remotely match anything assert f('boogie_monster') == 'c_boogie_monster' + + +@pytest.fixture +def fake_mixed_elastic(): + fake_tsl = openmc.data.ThermalScattering("c_D_in_7LiD", 1.9968, 4.9, [0.0253]) + fake_tsl.nuclides = ['H2'] + + # Create elastic reaction + bragg_edges = [0.00370672, 0.00494229, 0.00988458, 0.01359131, 0.01482688, + 0.01976918, 0.02347589, 0.02471147, 0.02965376, 0.03336048, + 0.03953834, 0.04324506, 0.04448063, 0.04942292, 0.05312964, + 0.05436522, 0.05930751, 0.06301423, 0.0642498 , 0.06919209, + 0.07289881, 0.07907667, 0.08278339, 0.08401896, 0.08896126, + 0.09266798, 0.09390355, 0.09884584, 0.1025526 , 0.1037882 , + 0.1087305 , 0.1124372 , 0.1186151 , 0.1223218 , 0.1235574 , + 0.1284997 , 0.1322064 , 0.133442 , 0.142091 , 0.1433266 , + 0.1482688 , 0.1519756 , 0.1581534 , 0.1618601 , 0.1630957 , + 0.168038 , 0.1717447 , 0.1729803 , 0.1779226 , 0.1816293 , + 0.1828649 , 0.1878072 , 0.1915139 , 0.1976918 , 0.2026341 , + 0.2075763 , 0.2125186 , 0.2174609 , 0.2224032 , 0.2273455 , + 0.2421724 , 0.2471147 , 0.252057 , 0.2569993 , 0.2619415 , + 0.2668838 , 0.2767684 , 0.2817107 , 0.2915953 , 0.3064222 , + 0.3261913 , 0.366965] + factors = [0.00375735, 0.01386287, 0.02595574, 0.02992438, 0.03549502, + 0.03855745, 0.04058831, 0.04986305, 0.05703106, 0.05855471, + 0.06078031, 0.06212291, 0.06656602, 0.06930339, 0.0697072 , + 0.07201456, 0.07263853, 0.07313129, 0.07465531, 0.07714482, + 0.07759976, 0.077809 , 0.07790282, 0.07927957, 0.08013058, + 0.08026637, 0.08073475, 0.08112202, 0.08123039, 0.08187171, + 0.08213756, 0.08218236, 0.08236572, 0.08240729, 0.08259795, + 0.08297893, 0.08300455, 0.08314566, 0.08315611, 0.08337715, + 0.08350026, 0.08350663, 0.08352815, 0.08353776, 0.0836098 , + 0.08367017, 0.08367361, 0.0837242 , 0.08375069, 0.08375227, + 0.08377006, 0.08381488, 0.08381644, 0.08382698, 0.08386266, + 0.08387756, 0.08388445, 0.08388974, 0.08390341, 0.08391088, + 0.08391695, 0.08392361, 0.08392684, 0.08392818, 0.08393161, + 0.08393546, 0.08393685, 0.08393801, 0.08393976, 0.08394167, + 0.08394288, 0.08394398] + coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors) + incoherent_xs = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672]) + elastic_xs = {'294K': openmc.data.Sum((coherent_xs, incoherent_xs))} + coherent_dist = openmc.data.CoherentElasticAE(coherent_xs) + incoherent_dist = openmc.data.IncoherentElasticAEDiscrete([ + [-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6] + ]) + elastic_dist = {'294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist)} + fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) + + # Create inelastic reaction + inelastic_xs = {'294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35])} + breakpoints = [3] + interpolation = [2] + energy = [1.0e-5, 4.3e-2, 4.9] + energy_out = [ + openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]), + ] + for eout in energy_out: + eout.normalize() + eout.c = eout.cdf() + discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8) + discrete.c = discrete.cdf()[1:] + mu = [[discrete]*4]*3 + inelastic_dist = {'294K': openmc.data.IncoherentInelasticAE( + breakpoints, interpolation, energy, energy_out, mu)} + inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist) + fake_tsl.inelastic = inelastic + + return fake_tsl + + +def test_mixed_elastic(fake_mixed_elastic, run_in_tmpdir): + # Write data to HDF5 and then read back + original = fake_mixed_elastic + original.export_to_hdf5('c_D_in_7LiD.h5') + copy = openmc.data.ThermalScattering.from_hdf5('c_D_in_7LiD.h5') + + # Make sure data did not change as a result of HDF5 writing/reading + assert original == copy + + # Create modified cross_sections.xml file that includes the above data + xs = openmc.data.DataLibrary.from_xml() + xs.register_file('c_D_in_7LiD.h5') + xs.export_to_xml('cross_sections_mixed.xml') + + # Create a minimal model that includes the new data and run it + mat = openmc.Material() + mat.add_nuclide('H2', 1.0) + mat.add_nuclide('Li7', 1.0) + mat.set_density('g/cm3', 1.0) + mat.add_s_alpha_beta('c_D_in_7LiD') + sph = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.materials = openmc.Materials([mat]) + model.materials.cross_sections = "cross_sections_mixed.xml" + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.Source( + energy=openmc.stats.Discrete([3.0], [1.0]) # 3 eV source + ) + model.run() diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index a0f18f851..b57242983 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -807,3 +807,46 @@ 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) + assert len(particles) == 10 + 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) + assert len(other_particles) == 10 + 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() diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index fd0e4f624..2c5cc2874 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -25,6 +25,52 @@ def test_add_nuclide(): with pytest.raises(ValueError): m.add_nuclide('H1', 1.0, 'oa') +def test_add_components(): + """Test adding multipe elements or nuclides at once""" + m = openmc.Material() + components = {'H1': 2.0, + 'O16': 1.0, + 'Zr': 1.0, + 'O': 1.0, + 'U': {'percent': 1.0, + 'enrichment': 4.5}, + 'Li': {'percent': 1.0, + 'enrichment': 60.0, + 'enrichment_target': 'Li7'}, + 'H': {'percent': 1.0, + 'enrichment': 50.0, + 'enrichment_target': 'H2', + 'enrichment_type': 'wo'}} + m.add_components(components) + with pytest.raises(ValueError): + m.add_components({'U': {'percent': 1.0, + 'enrichment': 100.0}}) + with pytest.raises(ValueError): + m.add_components({'Pu': {'percent': 1.0, + 'enrichment': 3.0}}) + with pytest.raises(ValueError): + m.add_components({'U': {'percent': 1.0, + 'enrichment': 70.0, + 'enrichment_target':'U235'}}) + with pytest.raises(ValueError): + m.add_components({'He': {'percent': 1.0, + 'enrichment': 17.0, + 'enrichment_target': 'He6'}}) + with pytest.raises(ValueError): + m.add_components({'li': 1.0}) # should fail as 1st char is lowercase + with pytest.raises(ValueError): + m.add_components({'LI': 1.0}) # should fail as 2nd char is uppercase + with pytest.raises(ValueError): + m.add_components({'Xx': 1.0}) # should fail as Xx is not an element + with pytest.raises(ValueError): + m.add_components({'n': 1.0}) # check to avoid n for neutron being accepted + with pytest.raises(TypeError): + m.add_components({'H1': '1.0'}) + with pytest.raises(TypeError): + m.add_components({1.0: 'H1'}, percent_type = 'wo') + with pytest.raises(ValueError): + m.add_components({'H1': 1.0}, percent_type = 'oa') + def test_remove_nuclide(): """Test removing nuclides.""" @@ -49,7 +95,7 @@ def test_remove_elements(): assert m.nuclides[0].percent == 1.0 -def test_elements(): +def test_add_element(): """Test adding elements.""" m = openmc.Material() m.add_element('Zr', 1.0) @@ -74,7 +120,6 @@ def test_elements(): with pytest.raises(ValueError): m.add_element('n', 1.0) # check to avoid n for neutron being accepted - def test_elements_by_name(): """Test adding elements by name""" m = openmc.Material() @@ -441,7 +486,7 @@ def test_activity_of_tritium(): m1.add_nuclide("H3", 1) m1.set_density('g/cm3', 1) m1.volume = 1 - assert pytest.approx(m1.activity) == 3.559778e14 + assert pytest.approx(m1.activity) == 3.559778e14 def test_activity_of_metastable(): diff --git a/tests/unit_tests/test_region.py b/tests/unit_tests/test_region.py index 8537e3b80..086ce6401 100644 --- a/tests/unit_tests/test_region.py +++ b/tests/unit_tests/test_region.py @@ -208,3 +208,17 @@ def test_from_expression(reset): # Opening parenthesis immediately after halfspace r = openmc.Region.from_expression('1(2|-3)', surfs) assert str(r) == '(1 (2 | -3))' + + +def test_translate_inplace(): + sph = openmc.Sphere() + x = openmc.XPlane() + region = -sph & +x + + # Translating a region should produce new surfaces + region2 = region.translate((0.5, -6.7, 3.9), inplace=False) + assert str(region) != str(region2) + + # Translating a region in-place should *not* produce new surfaces + region3 = region.translate((0.5, -6.7, 3.9), inplace=True) + assert str(region) == str(region3)