diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1673af531..e65e630b1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,7 +20,6 @@ env: COVERALLS_PARALLEL: true GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - jobs: main: runs-on: ubuntu-20.04 @@ -61,10 +60,10 @@ jobs: python-version: 3.8 omp: n mpi: y - name: 'Python ${{ matrix.python-version }} (omp=${{ matrix.omp }}, + name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }}, mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }}, libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }} - vectfit=${{ matrix.vectfit }})' + vectfit=${{ matrix.vectfit }})" env: MPI: ${{ matrix.mpi }} @@ -78,24 +77,23 @@ jobs: steps: - uses: actions/checkout@v2 - - - name: Set up Python ${{ matrix.python-version }} + - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} - - - name: Environment Variables + + - name: Environment Variables run: | echo "DAGMC_ROOT=$HOME/DAGMC" echo "OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml" >> $GITHUB_ENV echo "OPENMC_ENDF_DATA=$HOME/endf-b-vii.1" >> $GITHUB_ENV - - - name: Apt dependencies + - name: Apt dependencies shell: bash run: | sudo apt -y update - sudo apt install -y libmpich-dev \ + sudo apt install -y libpng-dev \ + libmpich-dev \ libnetcdf-dev \ libpnetcdf-dev \ libhdf5-serial-dev \ @@ -105,23 +103,21 @@ jobs: sudo update-alternatives --set mpirun /usr/bin/mpirun.mpich sudo update-alternatives --set mpi-x86_64-linux-gnu /usr/include/x86_64-linux-gnu/mpich - - - name: install + - name: install shell: bash run: | echo "$HOME/NJOY2016/build" >> $GITHUB_PATH $GITHUB_WORKSPACE/tools/ci/gha-install.sh - - - name: before + + - name: before shell: bash run: $GITHUB_WORKSPACE/tools/ci/gha-before-script.sh - - - name: test + + - name: test shell: bash run: $GITHUB_WORKSPACE/tools/ci/gha-script.sh - - - name: after_success + - name: after_success shell: bash run: | cpp-coveralls -i src -i include --exclude-pattern "/usr/*" --dump cpp_cov.json @@ -131,8 +127,8 @@ jobs: needs: main runs-on: ubuntu-latest steps: - - name: Coveralls Finished - uses: coverallsapp/github-action@master - with: - github-token: ${{ secrets.github_token }} - parallel-finished: true + - name: Coveralls Finished + uses: coverallsapp/github-action@master + with: + github-token: ${{ secrets.github_token }} + parallel-finished: true diff --git a/.gitignore b/.gitignore index f328d7b9a..88e59895d 100644 --- a/.gitignore +++ b/.gitignore @@ -64,6 +64,7 @@ scripts/*.tar.* scripts/G4EMLOW*/ # Images +*.png *.ppm *.voxel *.vti @@ -103,4 +104,4 @@ CMakeSettings.json .vscode/ # Python pickle files -*.pkl \ No newline at end of file +*.pkl diff --git a/CMakeLists.txt b/CMakeLists.txt index 7326e7f35..d0c2569c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -79,6 +79,11 @@ if(libmesh) find_package(LIBMESH REQUIRED) endif() +#=============================================================================== +# libpng +#=============================================================================== +find_package(PNG) + #=============================================================================== # HDF5 for binary output #=============================================================================== @@ -214,8 +219,6 @@ add_subdirectory(vendor/xtensor) # GSL header-only library #=============================================================================== -set(GSL_LITE_OPT_INSTALL_COMPAT_HEADER ON CACHE BOOL - "Install MS-GSL compatibility header ") add_subdirectory(vendor/gsl-lite) # Make sure contract violations throw exceptions @@ -425,6 +428,11 @@ if(libmesh) target_link_libraries(libopenmc PkgConfig::LIBMESH) endif() +if (PNG_FOUND) + target_compile_definitions(libopenmc PRIVATE USE_LIBPNG) + target_link_libraries(libopenmc PNG::PNG) +endif() + #=============================================================================== # openmc executable #=============================================================================== diff --git a/Dockerfile b/Dockerfile index fe87f9d5a..7803b3458 100644 --- a/Dockerfile +++ b/Dockerfile @@ -68,8 +68,8 @@ RUN apt-get update -y && \ apt-get install -y \ python3-pip python-is-python3 wget git build-essential cmake \ mpich libmpich-dev libhdf5-serial-dev libhdf5-mpich-dev \ - imagemagick && \ - apt-get autoremove -y + libpng-dev && \ + apt-get autoremove # Update system-provided pip RUN pip install --upgrade pip diff --git a/LICENSE b/LICENSE index 848a23845..3523dd370 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,5 @@ -Copyright (c) 2011-2021 Massachusetts Institute of Technology and OpenMC contributors +Copyright (c) 2011-2021 Massachusetts Institute of Technology, UChicago Argonne +LLC, and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index fa7ed7cc1..b12e95f23 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -16,6 +16,8 @@ if(@libmesh@) pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.6.0 IMPORTED_TARGET) endif() +find_package(PNG) + if(NOT TARGET OpenMC::libopenmc) include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") endif() diff --git a/docs/source/conf.py b/docs/source/conf.py index 52d872abf..a10494260 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -62,7 +62,7 @@ master_doc = 'index' # General information about the project. project = 'OpenMC' -copyright = '2011-2021, Massachusetts Institute of Technology and OpenMC contributors' +copyright = '2011-2021, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index 1d04bb8fb..b04edca2b 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -406,3 +406,14 @@ Each ```` element can have the following attributes or sub-eleme A required string indicating the file to be loaded representing the DAGMC universe. *Default*: None + + + .. note:: A geometry.xml file containing only a DAGMC model for a file named `dagmc.h5m` (no CSG) + looks as follows + + .. code-block:: xml + + + + + diff --git a/docs/source/io_formats/plots.rst b/docs/source/io_formats/plots.rst index 966d2b802..e6b75eafc 100644 --- a/docs/source/io_formats/plots.rst +++ b/docs/source/io_formats/plots.rst @@ -10,9 +10,9 @@ of the plots.xml is simply ```` and any number output plots can be defined with ```` sub-elements. Two plot types are currently implemented in openMC: -* ``slice`` 2D pixel plot along one of the major axes. Produces a PPM image +* ``slice`` 2D pixel plot along one of the major axes. Produces a PNG image file. -* ``voxel`` 3D voxel data dump. Produces a binary file containing voxel xyz +* ``voxel`` 3D voxel data dump. Produces an HDF5 file containing voxel xyz position and cell or material id. @@ -68,20 +68,14 @@ sub-elements: :type: Keyword for type of plot to be produced. Currently only "slice" and "voxel" plots are implemented. The "slice" plot type creates 2D pixel maps saved in - the PPM file format. PPM files can be displayed in most viewers (e.g. the - default Gnome viewer, IrfanView, etc.). The "voxel" plot type produces a - binary datafile containing voxel grid positioning and the cell or material - (specified by the ``color`` tag) at the center of each voxel. These - datafiles can be processed into VTK files using the :ref:`scripts_voxel` - script provided with OpenMC, and subsequently viewed with a 3D viewer such - as VISIT or Paraview. See the :ref:`io_voxel` for information about the - datafile structure. + the PNG file format. The "voxel" plot type produces a binary datafile + containing voxel grid positioning and the cell or material (specified by the + ``color`` tag) at the center of each voxel. Voxel plot files can be + processed into VTK files using the :ref:`scripts_voxel` script provided with + OpenMC and subsequently viewed with a 3D viewer such as VISIT or Paraview. + See the :ref:`io_voxel` for information about the datafile structure. - .. note:: Since the PPM format is saved without any kind of compression, - the resulting file sizes can be quite large. Saving the image in - the PNG format can often times reduce the file size by orders of - magnitude without any loss of image quality. Likewise, - high-resolution voxel files produced by OpenMC can be quite large, + .. note:: High-resolution voxel files produced by OpenMC can be quite large, but the equivalent VTK files will be significantly smaller. *Default*: "slice" @@ -94,11 +88,6 @@ attribute or sub-element: directions for "slice" and "voxel" plots, respectively. Should be two or three integers separated by spaces. - .. warning:: The ``pixels`` input determines the output file size. For the - PPM format, 10 million pixels will result in a file just under - 30 MB in size. A 10 million voxel binary file will be around - 40 MB. - .. warning:: If the aspect ratio defined in ``pixels`` does not match the aspect ratio defined in ``width`` the plot may appear stretched or squeezed. diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6ea8cfc56..c797212c1 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -619,16 +619,17 @@ variable and whose sub-elements/attributes are as follows: :type: The type of the distribution. Valid options are "uniform", "discrete", - "tabular", "maxwell", and "watt". The "uniform" option produces variates - sampled from a uniform distribution over a finite interval. The "discrete" - option produces random variates that can assume a finite number of values - (i.e., a distribution characterized by a probability mass function). The - "tabular" option produces random variates sampled from a tabulated + "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces + variates sampled from a uniform distribution over a finite interval. The + "discrete" option produces random variates that can assume a finite number + of values (i.e., a distribution characterized by a probability mass function). + The "tabular" option produces random variates sampled from a tabulated distribution where the density function is either a histogram or linearly-interpolated between tabulated points. The "watt" option produces random variates is sampled from a Watt fission spectrum (only used for energies). The "maxwell" option produce variates sampled from a Maxwell - fission spectrum (only used for energies). + fission spectrum (only used for energies). The "mixture" option produces samples + from univariate sub-distributions with given probabilities. *Default*: None @@ -651,12 +652,22 @@ variable and whose sub-elements/attributes are as follows: .. note:: The above format should be used even when using the multi-group :ref:`energy_mode`. + :interpolation: For a "tabular" distribution, ``interpolation`` can be set to "histogram" or "linear-linear" thereby specifying how tabular points are to be interpolated. *Default*: histogram +:pair: + For a "mixture" distribution, this element provides a distribution and its corresponding probability. + + :probability: + An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution. + + :dist: + This sub-element of a ``pair`` element provides information on the corresponding univariate distribution. + ------------------------- ```` Element ------------------------- diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index ad7227369..7f95f2985 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -318,6 +318,10 @@ attributes/sub-elements: :dimension: The number of mesh cells in each direction. (For regular mesh only.) + :length_multiplier: + A multiplicative factor to apply to the mesh coordinates in all directions. + (For unstructured mesh only.) + :lower_left: The lower-left corner of the structured mesh. If only two coordinates are given, it is assumed that the mesh is an x-y mesh. (For regular mesh only.) diff --git a/docs/source/license.rst b/docs/source/license.rst index 40ab106ad..d8ea319aa 100644 --- a/docs/source/license.rst +++ b/docs/source/license.rst @@ -4,7 +4,8 @@ License Agreement ================= -Copyright © 2011-2021 Massachusetts Institute of Technology and OpenMC contributors +Copyright © 2011-2021 Massachusetts Institute of Technology, UChicago Argonne +LLC, and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/source/quickinstall.rst b/docs/source/quickinstall.rst index e416b3efd..e3138a53b 100644 --- a/docs/source/quickinstall.rst +++ b/docs/source/quickinstall.rst @@ -100,7 +100,7 @@ directly from the package manager: .. code-block:: sh - sudo apt install g++ cmake libhdf5-dev + sudo apt install g++ cmake libhdf5-dev libpng-dev After the packages have been installed, follow the instructions below for building and installing OpenMC from source. diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 5fc1fe639..3090ae69e 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -219,6 +219,15 @@ Prerequisites .. admonition:: Optional :class: note + * libpng_ official reference PNG library + + OpenMC's built-in plotting capabilities use the libpng library to produce + compressed PNG files. In the absence of this library, OpenMC will fallback + to writing PPM files, which are uncompressed and only supported by select + image viewers. libpng can be installed on Ddebian derivates with:: + + sudo apt install libpng-dev + * An MPI implementation for distributed-memory parallel runs To compile with support for parallel runs on a distributed-memory diff --git a/docs/source/usersguide/plots.rst b/docs/source/usersguide/plots.rst index 109a8b549..d57917a3d 100644 --- a/docs/source/usersguide/plots.rst +++ b/docs/source/usersguide/plots.rst @@ -80,26 +80,13 @@ assign them to a :class:`openmc.Plots` collection and export it to XML:: plots += [plot2, plot3] plots.export_to_xml() -To actually generate the plots, run the :func:`openmc.plot_geometry` -function. Alternatively, run the :ref:`scripts_openmc` executable with the -``--plot`` command-line flag. When that has finished, you will have one or more -``.ppm`` files, i.e., `portable pixmap -`_ files. On some Linux -distributions, these ``.ppm`` files are natively viewable. If you find that -you're unable to open them on your system (or you don't like the fact that they -are not compressed), you may want to consider converting them to another format. -This is easily accomplished with the ``convert`` command available on most Linux -distributions as part of the `ImageMagick -`_ package. (On Debian -derivatives: ``sudo apt install imagemagick``). Images are then converted like: - -.. code-block:: sh - - convert myplot.ppm myplot.png - -Alternatively, if you're working within a `Jupyter `_ -Notebook or QtConsole, you can use the :func:`openmc.plot_inline` to run OpenMC -in plotting mode and display the resulting plot within the notebook. +To actually generate the plots, run the :func:`openmc.plot_geometry` function. +Alternatively, run the :ref:`scripts_openmc` executable with the ``--plot`` +command-line flag. When that has finished, you will have one or more ``.png`` +files. Alternatively, if you're working within a `Jupyter +`_ Notebook or QtConsole, you can use the +:func:`openmc.plot_inline` to run OpenMC in plotting mode and display the +resulting plot within the notebook. .. _usersguide_voxel: diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 5f6d33c97..d56fd5f81 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -10,7 +10,7 @@ #include "hdf5.h" #include "pugixml.hpp" -#include +#include #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr @@ -244,14 +244,14 @@ public: explicit CSGCell(pugi::xml_node cell_node); - bool contains(Position r, Direction u, int32_t on_surface) const; + bool contains(Position r, Direction u, int32_t on_surface) const override; std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const; + Position r, Direction u, int32_t on_surface, Particle* p) const override; void to_hdf5_inner(hid_t group_id) const override; - BoundingBox bounding_box() const; + BoundingBox bounding_box() const override; protected: bool contains_simple(Position r, Direction u, int32_t on_surface) const; diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 44fd985d2..8cdd78f78 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -24,6 +24,14 @@ public: virtual double sample(uint64_t* seed) const = 0; }; +using UPtrDist = unique_ptr; + +//! Return univariate probability distribution specified in XML file +//! \param[in] node XML node representing distribution +//! \return Unique pointer to distribution +UPtrDist distribution_from_xml(pugi::xml_node node); + + //============================================================================== //! A discrete distribution (probability mass function) //============================================================================== @@ -222,12 +230,27 @@ private: vector x_; //! Possible outcomes }; -using UPtrDist = unique_ptr; +//============================================================================== +//! Mixture distribution +//============================================================================== + +class Mixture : public Distribution { +public: + explicit Mixture(pugi::xml_node node); + + //! Sample a value from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled value + double sample(uint64_t* seed) const; + +private: + // Storrage for probability + distribution + using DistPair = std::pair; + + vector + distribution_; //!< sub-distributions + cummulative probabilities +}; -//! Return univariate probability distribution specified in XML file -//! \param[in] node XML node representing distribution -//! \return Unique pointer to distribution -UPtrDist distribution_from_xml(pugi::xml_node node); } // namespace openmc diff --git a/include/openmc/material.h b/include/openmc/material.h index 5fc52eba8..709d20573 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -6,7 +6,7 @@ #include "pugixml.hpp" #include "xtensor/xtensor.hpp" -#include +#include #include #include "openmc/bremsstrahlung.h" diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index d3ca900cf..d34ad2087 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -331,6 +331,15 @@ public: true}; //!< Write tallies onto the unstructured mesh at the end of a run std::string filename_; //!< Path to unstructured mesh file +protected: + //! Set the length multiplier to apply to each point in the mesh + void set_length_multiplier(const double length_multiplier); + + // Data members + double length_multiplier_ { + 1.0}; //!< Constant multiplication factor to apply to mesh coordinates + bool specified_length_multiplier_ {false}; + private: //! Setup method for the mesh. Builds data structures, //! sets up element mapping, creates bounding boxes, etc. @@ -344,7 +353,7 @@ public: // Constructors MOABMesh() = default; MOABMesh(pugi::xml_node); - MOABMesh(const std::string& filename); + MOABMesh(const std::string& filename, double length_multiplier = 1.0); MOABMesh(std::shared_ptr external_mbi); // Overridden Methods @@ -492,7 +501,7 @@ class LibMesh : public UnstructuredMesh { public: // Constructors LibMesh(pugi::xml_node node); - LibMesh(const std::string& filename); + LibMesh(const std::string& filename, double length_multiplier = 1.0); // Overridden Methods void bins_crossed(Position r0, Position r1, const Direction& u, diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index f961e9808..b05c76bbb 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -7,7 +7,7 @@ #include #include // for pair -#include +#include #include #include "openmc/array.h" diff --git a/include/openmc/photon.h b/include/openmc/photon.h index bf64fa50c..09783fdb4 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -7,7 +7,7 @@ #include "openmc/vector.h" #include "xtensor/xtensor.hpp" -#include +#include #include #include diff --git a/include/openmc/plot.h b/include/openmc/plot.h index aa15277aa..650b7e16a 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -30,9 +30,7 @@ namespace model { extern std::unordered_map plot_map; //!< map of plot ids to index extern vector plots; //!< Plot instance container -extern uint64_t - plotter_prn_seeds[N_STREAMS]; // Random number seeds used for plotter -extern int plotter_stream; // Stream index used by the plotter +extern uint64_t plotter_seed; // Stream index used by the plotter } // namespace model @@ -239,11 +237,18 @@ public: //! \param[out] image data associated with the plot object void draw_mesh_lines(Plot const& pl, ImageData& data); -//! Write a ppm image to file using a plot object's image data +//! Write a PPM image using a plot object's image data //! \param[in] plot object //! \param[out] image data associated with the plot object void output_ppm(Plot const& pl, const ImageData& data); +#ifdef USE_LIBPNG +//! Write a PNG image using a plot object's image data +//! \param[in] plot object +//! \param[out] image data associated with the plot object +void output_png(Plot const& pl, const ImageData& data); +#endif + //! Initialize a voxel file //! \param[in] id of an open hdf5 file //! \param[in] dimensions of the voxel file (dx, dy, dz) @@ -274,9 +279,12 @@ void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace); //! Read plot specifications from a plots.xml file void read_plots_xml(); -//! Create a ppm image for a plot object +//! Clear memory +void free_memory_plot(); + +//! Create an image for a plot object //! \param[in] plot object -void create_ppm(Plot const& pl); +void create_image(Plot const& pl); //! Create an hdf5 voxel file for a plot object //! \param[in] plot object diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 6b276db20..46705acf5 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -7,7 +7,7 @@ #include #include "hdf5.h" -#include +#include #include "openmc/reaction_product.h" #include "openmc/vector.h" diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 2d0261960..f14f6c15f 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -32,7 +32,8 @@ extern "C" bool cmfd_run; //!< is a CMFD run? extern bool delayed_photon_scaling; //!< Scale fission photon yield to include delayed extern "C" bool entropy_on; //!< calculate Shannon entropy? -extern bool event_based; //!< use event-based mode (instead of history-based) +extern "C" bool + event_based; //!< use event-based mode (instead of history-based) extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular? extern bool material_cell_offsets; //!< create material cells offsets? extern "C" bool output_summary; //!< write summary.h5? diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index b0145f0e1..2eecfe0e4 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -6,7 +6,7 @@ #include #include "pugixml.hpp" -#include +#include #include "openmc/constants.h" #include "openmc/hdf5_interface.h" diff --git a/include/openmc/tallies/filter_azimuthal.h b/include/openmc/tallies/filter_azimuthal.h index 4cda9163f..2272b500a 100644 --- a/include/openmc/tallies/filter_azimuthal.h +++ b/include/openmc/tallies/filter_azimuthal.h @@ -4,7 +4,7 @@ #include "openmc/vector.h" #include -#include +#include #include "openmc/tallies/filter.h" diff --git a/include/openmc/tallies/filter_cell.h b/include/openmc/tallies/filter_cell.h index 6e2463f24..46d89811d 100644 --- a/include/openmc/tallies/filter_cell.h +++ b/include/openmc/tallies/filter_cell.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_cell_instance.h b/include/openmc/tallies/filter_cell_instance.h index d74850df7..4de3fcd29 100644 --- a/include/openmc/tallies/filter_cell_instance.h +++ b/include/openmc/tallies/filter_cell_instance.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include "openmc/cell.h" #include "openmc/tallies/filter.h" diff --git a/include/openmc/tallies/filter_collision.h b/include/openmc/tallies/filter_collision.h index 46c64c553..3724b06cf 100644 --- a/include/openmc/tallies/filter_collision.h +++ b/include/openmc/tallies/filter_collision.h @@ -1,7 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_COLLISIONS_H #define OPENMC_TALLIES_FILTER_COLLISIONS_H -#include +#include #include #include "openmc/tallies/filter.h" diff --git a/include/openmc/tallies/filter_delayedgroup.h b/include/openmc/tallies/filter_delayedgroup.h index 868807f05..72ffa1db5 100644 --- a/include/openmc/tallies/filter_delayedgroup.h +++ b/include/openmc/tallies/filter_delayedgroup.h @@ -1,7 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_DELAYEDGROUP_H #define OPENMC_TALLIES_FILTER_DELAYEDGROUP_H -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h index c820f28fb..000aaa28b 100644 --- a/include/openmc/tallies/filter_energy.h +++ b/include/openmc/tallies/filter_energy.h @@ -1,7 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_ENERGY_H #define OPENMC_TALLIES_FILTER_ENERGY_H -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_material.h b/include/openmc/tallies/filter_material.h index 8a67c1a75..f58fc9938 100644 --- a/include/openmc/tallies/filter_material.h +++ b/include/openmc/tallies/filter_material.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h index de934d156..5299f6dd4 100644 --- a/include/openmc/tallies/filter_mu.h +++ b/include/openmc/tallies/filter_mu.h @@ -1,7 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_MU_H #define OPENMC_TALLIES_FILTER_MU_H -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_polar.h b/include/openmc/tallies/filter_polar.h index 23c00e619..e06aca1e0 100644 --- a/include/openmc/tallies/filter_polar.h +++ b/include/openmc/tallies/filter_polar.h @@ -3,7 +3,7 @@ #include -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_sph_harm.h b/include/openmc/tallies/filter_sph_harm.h index 498590d3c..5f5bf84f2 100644 --- a/include/openmc/tallies/filter_sph_harm.h +++ b/include/openmc/tallies/filter_sph_harm.h @@ -3,7 +3,7 @@ #include -#include +#include #include "openmc/tallies/filter.h" diff --git a/include/openmc/tallies/filter_surface.h b/include/openmc/tallies/filter_surface.h index 368fd09d0..358963fde 100644 --- a/include/openmc/tallies/filter_surface.h +++ b/include/openmc/tallies/filter_surface.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/filter_universe.h b/include/openmc/tallies/filter_universe.h index afe583b2d..fde0b6397 100644 --- a/include/openmc/tallies/filter_universe.h +++ b/include/openmc/tallies/filter_universe.h @@ -4,7 +4,7 @@ #include #include -#include +#include #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 39e0b0af2..13b8317ee 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -10,7 +10,7 @@ #include "pugixml.hpp" #include "xtensor/xfixed.hpp" #include "xtensor/xtensor.hpp" -#include +#include #include #include diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index 33cc7343a..db96f250f 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -9,7 +9,7 @@ #include "pugixml.hpp" #include "xtensor/xtensor.hpp" -#include +#include #include namespace openmc { diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 8a93af514..f99389069 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -57,8 +57,8 @@ Indicates the default path to an HDF5 file that contains multi-group cross section libraries if the user has not specified the tag in .I materials.xml\fP. .SH LICENSE -Copyright \(co 2011-2021 Massachusetts Institute of Technology and OpenMC -contributors. +Copyright \(co 2011-2021 Massachusetts Institute of Technology, UChicago +Argonne LLC, and OpenMC contributors. .PP Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/openmc/__init__.py b/openmc/__init__.py index 8985099ba..d118e6a05 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -34,4 +34,5 @@ from . import examples # Import a few names from the model module from openmc.model import rectangular_prism, hexagonal_prism, Model + __version__ = '0.13.0-dev' diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 3e7d9599c..a0402b47d 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -57,8 +57,13 @@ def get_decay_modes(value): List of successive decays, e.g. ('beta-', 'neutron') """ - return [_DECAY_MODES[int(x)][0] for x in - str(value).strip('0').replace('.', '')] + if int(value) == 10: + # The logic below would treat 10.0 as [1, 0] rather than [10] as it + # should, so we handle this case separately + return ['unknown'] + else: + return [_DECAY_MODES[int(x)][0] for x in + str(value).strip('0').replace('.', '')] class FissionProductYields(EqualityMixin): diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 0ef77695b..305edaf4f 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -5,95 +5,74 @@ import shutil from subprocess import Popen, PIPE, STDOUT, CalledProcessError import tempfile from pathlib import Path +import warnings from . import endf +import openmc.data -# For a given MAT number, give a name for the ACE table and a list of ZAID -# identifiers. This is based on Appendix C in the ENDF manual. +# For a given material, give a name for the ACE table and a list of ZAID +# identifiers. ThermalTuple = namedtuple('ThermalTuple', ['name', 'zaids', 'nmix']) _THERMAL_DATA = { - 1: ThermalTuple('hh2o', [1001], 1), - 2: ThermalTuple('parah', [1001], 1), - 3: ThermalTuple('orthoh', [1001], 1), - 5: ThermalTuple('hyh2', [1001], 1), - 7: ThermalTuple('hzrh', [1001], 1), - 8: ThermalTuple('hcah2', [1001], 1), - 10: ThermalTuple('hice', [1001], 1), - 11: ThermalTuple('dd2o', [1002], 1), - 12: ThermalTuple('parad', [1002], 1), - 13: ThermalTuple('orthod', [1002], 1), - 14: ThermalTuple('dice', [1002], 1), - 26: ThermalTuple('be', [4009], 1), - 27: ThermalTuple('bebeo', [4009], 1), - 28: ThermalTuple('bebe2c', [4009], 1), - 30: ThermalTuple('graph', [6000, 6012, 6013], 1), - 31: ThermalTuple('grph10', [6000, 6012, 6013], 1), - 32: ThermalTuple('grph30', [6000, 6012, 6013], 1), - 33: ThermalTuple('lch4', [1001], 1), - 34: ThermalTuple('sch4', [1001], 1), - 35: ThermalTuple('sch4p2', [1001], 1), - 37: ThermalTuple('hch2', [1001], 1), - 38: ThermalTuple('mesi00', [1001], 1), - 39: ThermalTuple('lucite', [1001], 1), - 40: ThermalTuple('benz', [1001, 6000, 6012], 2), - 42: ThermalTuple('tol00', [1001], 1), - 43: ThermalTuple('sisic', [14028, 14029, 14030], 1), - 44: ThermalTuple('csic', [6000, 6012, 6013], 1), - 45: ThermalTuple('ouo2', [8016, 8017, 8018], 1), - 46: ThermalTuple('obeo', [8016, 8017, 8018], 1), - 47: ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3), - 48: ThermalTuple('osap00', [92238], 1), - 49: ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3), - 50: ThermalTuple('oice', [8016, 8017, 8018], 1), - 51: ThermalTuple('od2o', [8016, 8017, 8018], 1), - 52: ThermalTuple('mg24', [12024], 1), - 53: ThermalTuple('al27', [13027], 1), - 55: ThermalTuple('yyh2', [39089], 1), - 56: ThermalTuple('fe56', [26056], 1), - 58: ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1), - 59: ThermalTuple('si00', [14028], 1), - 60: ThermalTuple('asap00', [13027], 1), - 71: ThermalTuple('n-un', [7014, 7015], 1), - 72: ThermalTuple('u-un', [92238], 1), - 75: ThermalTuple('uuo2', [8016, 8017, 8018], 1), + 'c_Al27': ThermalTuple('al27', [13027], 1), + 'c_Al_in_Al2O3': ThermalTuple('asap00', [13027], 1), + 'c_Be': ThermalTuple('be', [4009], 1), + 'c_Be_in_BeO': ThermalTuple('bebeo', [4009], 1), + 'c_Be_in_Be2C': ThermalTuple('bebe2c', [4009], 1), + 'c_Be_in_FLiBe': ThermalTuple('beflib', [4009], 1), + 'c_C6H6': ThermalTuple('benz', [1001, 6000, 6012], 2), + 'c_C_in_SiC': ThermalTuple('csic', [6000, 6012, 6013], 1), + 'c_Ca_in_CaH2': ThermalTuple('cacah2', [20040, 20042, 20043, 20044, 20046, 20048], 1), + 'c_D_in_D2O': ThermalTuple('dd2o', [1002], 1), + 'c_D_in_D2O_solid': ThermalTuple('dice', [1002], 1), + 'c_F_in_FLiBe': ThermalTuple('fflibe', [9019], 1), + 'c_Fe56': ThermalTuple('fe56', [26056], 1), + 'c_Graphite': ThermalTuple('graph', [6000, 6012, 6013], 1), + 'c_Graphite_10p': ThermalTuple('grph10', [6000, 6012, 6013], 1), + 'c_Graphite_30p': ThermalTuple('grph30', [6000, 6012, 6013], 1), + 'c_H_in_C5O2H8': ThermalTuple('lucite', [1001], 1), + 'c_H_in_CaH2': ThermalTuple('hcah2', [1001], 1), + 'c_H_in_CH2': ThermalTuple('hch2', [1001], 1), + 'c_H_in_CH4_liquid': ThermalTuple('lch4', [1001], 1), + 'c_H_in_CH4_solid': ThermalTuple('sch4', [1001], 1), + 'c_H_in_CH4_solid_phase_II': ThermalTuple('sch4p2', [1001], 1), + 'c_H_in_H2O': ThermalTuple('hh2o', [1001], 1), + 'c_H_in_H2O_solid': ThermalTuple('hice', [1001], 1), + 'c_H_in_HF': ThermalTuple('hhf', [1001], 1), + 'c_H_in_Mesitylene': ThermalTuple('mesi00', [1001], 1), + 'c_H_in_ParaffinicOil': ThermalTuple('hparaf', [1001], 1), + 'c_H_in_Toluene': ThermalTuple('tol00', [1001], 1), + 'c_H_in_UH3': ThermalTuple('huh3', [1001], 1), + 'c_H_in_YH2': ThermalTuple('hyh2', [1001], 1), + 'c_H_in_ZrH': ThermalTuple('hzrh', [1001], 1), + 'c_H_in_ZrH2': ThermalTuple('hzrh2', [1001], 1), + 'c_H_in_ZrHx': ThermalTuple('hzrhx', [1001], 1), + 'c_Li_in_FLiBe': ThermalTuple('liflib', [3006, 3007], 1), + 'c_Mg24': ThermalTuple('mg24', [12024], 1), + 'c_N_in_UN': ThermalTuple('n-un', [7014, 7015], 1), + 'c_O_in_Al2O3': ThermalTuple('osap00', [92238], 1), + 'c_O_in_BeO': ThermalTuple('obeo', [8016, 8017, 8018], 1), + 'c_O_in_D2O': ThermalTuple('od2o', [8016, 8017, 8018], 1), + 'c_O_in_H2O_solid': ThermalTuple('oice', [8016, 8017, 8018], 1), + 'c_O_in_UO2': ThermalTuple('ouo2', [8016, 8017, 8018], 1), + 'c_ortho_D': ThermalTuple('orthod', [1002], 1), + 'c_ortho_H': ThermalTuple('orthoh', [1001], 1), + 'c_para_D': ThermalTuple('parad', [1002], 1), + 'c_para_H': ThermalTuple('parah', [1001], 1), + 'c_Si28': ThermalTuple('si00', [14028], 1), + 'c_Si_in_SiC': ThermalTuple('sisic', [14028, 14029, 14030], 1), + 'c_SiO2_alpha': ThermalTuple('sio2-a', [8016, 8017, 8018, 14028, 14029, 14030], 3), + 'c_SiO2_beta': ThermalTuple('sio2-b', [8016, 8017, 8018, 14028, 14029, 14030], 3), + 'c_U_in_UN': ThermalTuple('u-un', [92238], 1), + 'c_U_in_UO2': ThermalTuple('uuo2', [8016, 8017, 8018], 1), + 'c_Y_in_YH2': ThermalTuple('yyh2', [39089], 1), + 'c_Zr_in_ZrH': ThermalTuple('zrzrh', [40000, 40090, 40091, 40092, 40094, 40096], 1), + 'c_Zr_in_ZrH2': ThermalTuple('zrzrh2', [40000, 40090, 40091, 40092, 40094, 40096], 1), + 'c_Zr_in_ZrHx': ThermalTuple('zrzrhx', [40000, 40090, 40091, 40092, 40094, 40096], 1), } -def _get_thermal_data(ev, mat): - """Return appropriate ThermalTuple, accounting for bugs.""" - - # JEFF assigns MAT=59 to Ca in CaH2 (which is supposed to be silicon). - if ev.info['library'][0] == 'JEFF': - if ev.material == 59: - if 'CaH2' in ''.join(ev.info['description']): - zaids = [20040, 20042, 20043, 20044, 20046, 20048] - return ThermalTuple('cacah2', zaids, 1) - - # Before ENDF/B-VIII.0, crystalline graphite was MAT=31 - if ev.info['library'] != ('ENDF/B', 8, 0): - if ev.material == 31: - return _THERMAL_DATA[30] - - # ENDF/B incorrectly assigns MAT numbers for UO2 - # - # Material | ENDF Manual | VII.0 | VII.1 | VIII.0 - # ---------|-------------|-------|-------|------- - # O in UO2 | 45 | 75 | 75 | 75 - # U in UO2 | 75 | 76 | 48 | 48 - if ev.info['library'][0] == 'ENDF/B': - if ev.material == 75: - return _THERMAL_DATA[45] - version = ev.info['library'][1:] - if version in ((7, 1), (8, 0)) and ev.material == 48: - return _THERMAL_DATA[75] - if version == (7, 0) and ev.material == 76: - return _THERMAL_DATA[75] - - # If not a problematic material, use the dictionary as is - return _THERMAL_DATA[mat] - - _TEMPLATE_RECONR = """ reconr / %%%%%%%%%%%%%%%%%%% Reconstruct XS for neutrons %%%%%%%%%%%%%%%%%%%%%%% {nendf} {npendf} @@ -170,7 +149,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% {nendf} {nthermal_acer_in} 0 {nace} {ndir} 2 0 1 .{ext}/ '{library}: {zsymam_thermal} processed by NJOY'/ -{mat} {temperature} '{data.name}' / +{mat} {temperature} '{data.name}' {nza} / {zaids} / 222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/ """ @@ -445,7 +424,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, def make_ace_thermal(filename, filename_thermal, temperatures=None, ace='ace', xsdir=None, output_dir=None, error=0.001, - iwt=2, evaluation=None, evaluation_thermal=None, **kwargs): + iwt=2, evaluation=None, evaluation_thermal=None, + table_name=None, zaids=None, nmix=None, **kwargs): """Generate thermal scattering ACE file from ENDF files Parameters @@ -476,6 +456,12 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None, evaluation_thermal : openmc.data.endf.Evaluation, optional If the ENDF thermal scattering sublibrary file contains multiple material evaluations, this argument indicates which evaluation to use. + table_name : str, optional + Name to assign to ACE table + zaids : list of int, optional + ZAIDs that the thermal scattering data applies to + nmix : int, optional + Number of atom types in mixed moderator **kwargs Keyword arguments passed to :func:`openmc.data.njoy.run` @@ -499,11 +485,23 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None, ev_thermal = (evaluation_thermal if evaluation_thermal is not None else endf.Evaluation(filename_thermal)) mat_thermal = ev_thermal.material - zsymam_thermal = ev_thermal.target['zsymam'] + zsymam_thermal = ev_thermal.target['zsymam'].strip() - # Determine name, isotopes based on MAT number - data = _get_thermal_data(ev_thermal, mat_thermal) - zaids = ' '.join(str(zaid) for zaid in data.zaids[:3]) + # Determine name, isotopes, and number of atom types + if table_name and zaids and nmix: + data = ThermalTuple(table_name, zaids, nmix) + else: + with warnings.catch_warnings(record=True) as w: + proper_name = openmc.data.get_thermal_name(zsymam_thermal) + if w: + raise RuntimeError( + f"Thermal scattering material {zsymam_thermal} not " + "recognized. Please contact OpenMC developers at " + "https://openmc.discourse.group.") + data = _THERMAL_DATA[proper_name] + + zaids = ' '.join(str(zaid) for zaid in data.zaids) + nza = len(data.zaids) # Determine name of library library = '{}-{}.{}'.format(*ev_thermal.info['library']) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index d0adbc4fe..b807a9f59 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -28,52 +28,62 @@ from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, _THERMAL_NAMES = { - 'c_Al27': ('al', 'al27', 'al-27'), - 'c_Al_in_Sapphire': ('asap00', 'asap'), - 'c_Be': ('be', 'be-metal', 'be-met', 'be00'), + 'c_Al27': ('al', 'al27', 'al-27', '13-al- 27'), + 'c_Al_in_Al2O3': ('asap00', 'asap', 'al(al2o3)'), + 'c_Be': ('be', 'be-metal', 'be-met', 'be00', 'be-metal', 'be metal', '4-be'), 'c_BeO': ('beo',), - 'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o', 'bbeo00'), + 'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o', 'bbeo00', 'be(beo)'), 'c_Be_in_Be2C': ('bebe2c',), - 'c_C6H6': ('benz', 'c6h6'), - 'c_C_in_SiC': ('csic', 'c-sic'), - 'c_Ca_in_CaH2': ('cah', 'cah00', 'cacah2'), - 'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw', 'dhw00'), + 'c_Be_in_FLiBe': ('beflib', 'be(flibe)'), + 'c_C6H6': ('benz', 'c6h6', 'benzine'), + 'c_C_in_SiC': ('csic', 'c-sic', 'c(3c-sic)'), + 'c_Ca_in_CaH2': ('cah', 'cah00', 'cacah2', 'ca(cah2)'), + 'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw', 'dhw00', 'd(d2o)'), 'c_D_in_D2O_ice': ('dice',), - 'c_Fe56': ('fe', 'fe56', 'fe-56'), - 'c_Graphite': ('graph', 'grph', 'gr', 'gr00'), - 'c_Graphite_10p': ('grph10',), - 'c_Graphite_30p': ('grph30',), - 'c_H_in_CaH2': ('hcah2', 'hca00'), - 'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00'), - 'c_H_in_CH4_liquid': ('lch4', 'lmeth'), - 'c_H_in_CH4_solid': ('sch4', 'smeth'), + 'c_F_in_FLiBe': ('fflibe', 'f(flibe)'), + 'c_Fe56': ('fe', 'fe56', 'fe-56', '26-fe- 56'), + 'c_Graphite': ('graph', 'grph', 'gr', 'gr00', 'graphite'), + 'c_Graphite_10p': ('grph10', '10p graphit'), + 'c_Graphite_30p': ('grph30', '30p graphit'), + 'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci', 'h(lucite)'), + 'c_H_in_CaH2': ('hcah2', 'hca00', 'h(cah2)'), + 'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00', 'h(ch2)'), + 'c_H_in_CH4_liquid': ('lch4', 'lmeth', 'l-ch4'), + 'c_H_in_CH4_solid': ('sch4', 'smeth', 's-ch4'), 'c_H_in_CH4_solid_phase_II': ('sch4p2',), - 'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw', 'lw00'), - 'c_H_in_H2O_solid': ('hice', 'h-ice', 'ice00'), - 'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci'), - 'c_H_in_Mesitylene': ('mesi00', 'mesi'), - 'c_H_in_Toluene': ('tol00', 'tol'), - 'c_H_in_YH2': ('hyh2', 'h-yh2'), - 'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr', 'hzr00'), - 'c_Mg24': ('mg', 'mg24', 'mg00'), - 'c_O_in_Sapphire': ('osap00', 'osap'), - 'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be', 'obeo00'), - 'c_O_in_D2O': ('od2o', 'o-d2o', 'ohw00'), - 'c_O_in_H2O_ice': ('oice', 'o-ice'), - 'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u', 'ouo200'), - 'c_N_in_UN': ('n-un',), - 'c_ortho_D': ('orthod', 'orthoD', 'dortho', 'od200', 'ortod'), - 'c_ortho_H': ('orthoh', 'orthoH', 'hortho', 'oh200', 'ortoh'), - 'c_Si28': ('si00', 'sili'), - 'c_Si_in_SiC': ('sisic', 'si-sic'), - 'c_SiO2_alpha': ('sio2', 'sio2a'), - 'c_SiO2_beta': ('sio2b',), - 'c_para_D': ('parad', 'paraD', 'dpara', 'pd200'), - 'c_para_H': ('parah', 'paraH', 'hpara', 'ph200'), - 'c_U_in_UN': ('u-un',), - 'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2', 'uuo200'), - 'c_Y_in_YH2': ('yyh2', 'y-yh2'), - 'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h') + 'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw', 'lw00', 'h(h2o)'), + 'c_H_in_H2O_solid': ('hice', 'h-ice', 'ice00', 'h(ice-ih)', 'h(ice)'), + 'c_H_in_HF': ('hhf', 'h(hf)'), + 'c_H_in_Mesitylene': ('mesi00', 'mesi', 'mesi-phii'), + 'c_H_in_ParaffinicOil': ('hparaf', 'h(paraffin'), + 'c_H_in_Toluene': ('tol00', 'tol', 'tolue-phii'), + 'c_H_in_UH3': ('huh3', 'h(uh3)'), + 'c_H_in_YH2': ('hyh2', 'h-yh2', 'h(yh2)'), + 'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr', 'hzr00', 'h(zrh)'), + 'c_H_in_ZrH2': ('hzrh2', 'h(zrh2)'), + 'c_H_in_ZrHx': ('hzrhx', 'h(zrhx)'), + 'c_Li_in_FLiBe': ('liflib', 'li(flibe)'), + 'c_Mg24': ('mg', 'mg24', 'mg00', '24-mg'), + 'c_N_in_UN': ('n-un', 'n(un)', 'n(un) l'), + 'c_O_in_Al2O3': ('osap00', 'osap', 'o(al2o3)'), + 'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be', 'obeo00', 'o(beo)'), + 'c_O_in_D2O': ('od2o', 'o-d2o', 'ohw00', 'o(d2o)'), + 'c_O_in_H2O_solid': ('oice', 'o-ice', 'o(ice-ih)'), + 'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u', 'ouo200', 'o(uo2)'), + 'c_ortho_D': ('orthod', 'orthoD', 'dortho', 'od200', 'ortod', 'ortho-d'), + 'c_ortho_H': ('orthoh', 'orthoH', 'hortho', 'oh200', 'ortoh', 'ortho-h'), + 'c_para_D': ('parad', 'paraD', 'dpara', 'pd200', 'para-d'), + 'c_para_H': ('parah', 'paraH', 'hpara', 'ph200', 'para-h'), + 'c_Si28': ('si00', 'sili', 'si'), + 'c_Si_in_SiC': ('sisic', 'si-sic', 'si(3c-sic)'), + 'c_SiO2_alpha': ('sio2', 'sio2a', 'sio2alpha'), + 'c_SiO2_beta': ('sio2b', 'sio2beta'), + 'c_U_in_UN': ('u-un', 'u(un)', 'u(un) l'), + 'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2', 'uuo200', 'u(uo2)'), + 'c_Y_in_YH2': ('yyh2', 'y-yh2', 'y(yh2)'), + 'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h', 'zr(zrh)'), + 'c_Zr_in_ZrH2': ('zrzrh2', 'zr(zrh2)'), + 'c_Zr_in_ZrHx': ('zrzrhx', 'zr(zrhx)'), } diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 00b29e319..b8c1cdfff 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -4,17 +4,6 @@ openmc.deplete A depletion front-end tool. """ -import sys -from unittest.mock import Mock - -from .dummy_comm import DummyCommunicator - -try: - from mpi4py import MPI - comm = MPI.COMM_WORLD -except ImportError: - MPI = Mock() - comm = DummyCommunicator() from .nuclide import * from .chain import * diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ca16926ea..b9198696d 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -22,7 +22,7 @@ from uncertainties import ufloat from openmc.data import DataLibrary from openmc.lib import MaterialFilter, Tally from openmc.checkvalue import check_type, check_greater_than -from . import comm +from openmc.mpi import comm from .results import Results from .chain import Chain from .results_list import ResultsList diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 6c3badc09..f22382923 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -10,7 +10,7 @@ import sys from numpy import dot, zeros, newaxis, asarray -from . import comm +from openmc.mpi import comm from openmc.checkvalue import check_type, check_greater_than from openmc.data import JOULE_PER_EV, REACTION_MT from openmc.lib import ( diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index fbf43a808..585deabba 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -581,3 +581,14 @@ class SILEQIIntegrator(SIIntegrator): proc_time += time1 + time2 return proc_time, [eos_conc, inter_conc], [res_bar] + + +integrator_by_name = { + 'cecm': CECMIntegrator, + 'predictor': PredictorIntegrator, + 'cf4': CF4Integrator, + 'epc_rk4': EPCRK4Integrator, + 'si_celi': SICELIIntegrator, + 'si_leqi': SILEQIIntegrator, + 'celi': CELIIntegrator, + 'leqi': LEQIIntegrator} diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index edba9748c..992248246 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -294,7 +294,8 @@ class Nuclide: for mode_type, daughter, br in self.decay_modes: mode_elem = ET.SubElement(elem, 'decay') mode_elem.set('type', mode_type) - mode_elem.set('target', daughter or "Nothing") + if daughter: + mode_elem.set('target', daughter) mode_elem.set('branching_ratio', str(br)) elem.set('reactions', str(len(self.reactions))) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 0e5cc9efa..29370eb8c 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -20,7 +20,7 @@ from uncertainties import ufloat import openmc from openmc.checkvalue import check_value import openmc.lib -from . import comm +from openmc.mpi import comm from .abc import TransportOperator, OperatorResult from .atom_number import AtomNumber from .reaction_rates import ReactionRates @@ -176,6 +176,9 @@ class Operator(TransportOperator): results are to be used. diff_burnable_mats : bool Whether to differentiate burnable materials with multiple instances + cleanup_when_done : bool + Whether to finalize and clear the shared library memory when the + depletion operation is complete. Defaults to clearing the library. """ _fission_helpers = { "average": AveragedFissionYieldHelper, @@ -202,6 +205,7 @@ class Operator(TransportOperator): self.settings = settings self.geometry = geometry self.diff_burnable_mats = diff_burnable_mats + self.cleanup_when_done = True # Reduce the chain before we create more materials if reduce_chain: @@ -317,6 +321,10 @@ class Operator(TransportOperator): # Reset results in OpenMC openmc.lib.reset() + # Update the number densities regardless of the source rate + self.number.set_density(vec) + self._update_materials() + # If the source rate is zero, return zero reaction rates without running # a transport solve if source_rate == 0.0: @@ -327,11 +335,7 @@ class Operator(TransportOperator): # Prevent OpenMC from complaining about re-creating tallies openmc.reset_auto_ids() - # Update status - self.number.set_density(vec) - - # Update material compositions and tally nuclides - self._update_materials() + # Update tally nuclides data in preparation for transport solve nuclides = self._get_tally_nuclides() self._rate_helper.nuclides = nuclides self._normalization_helper.nuclides = nuclides @@ -536,7 +540,8 @@ class Operator(TransportOperator): # Initialize OpenMC library comm.barrier() - openmc.lib.init(intracomm=comm) + if not openmc.lib.is_initialized: + openmc.lib.init(intracomm=comm) # Generate tallies in memory materials = [openmc.lib.materials[int(i)] @@ -554,7 +559,8 @@ class Operator(TransportOperator): def finalize(self): """Finalize a depletion simulation and release resources.""" - openmc.lib.finalize() + if self.cleanup_when_done: + openmc.lib.finalize() def _update_materials(self): """Updates material compositions in OpenMC on all processes.""" diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index ce8331e82..111c8e638 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -9,7 +9,7 @@ import copy import h5py import numpy as np -from . import comm, MPI +from openmc.mpi import comm, MPI from .reaction_rates import ReactionRates VERSION_RESULTS = (1, 1) diff --git a/openmc/deplete/dummy_comm.py b/openmc/dummy_comm.py similarity index 100% rename from openmc/deplete/dummy_comm.py rename to openmc/dummy_comm.py diff --git a/openmc/executor.py b/openmc/executor.py index e8f0de2c7..63653f4d0 100644 --- a/openmc/executor.py +++ b/openmc/executor.py @@ -3,6 +3,84 @@ from numbers import Integral import subprocess import openmc +from .plots import _get_plot_image + + +def _process_CLI_arguments(volume=False, geometry_debug=False, particles=None, + plot=False, restart_file=None, threads=None, + tracks=False, event_based=None, + openmc_exec='openmc', mpi_args=None): + """Converts user-readable flags in to command-line arguments to be run with + the OpenMC executable via subprocess. + + Parameters + ---------- + volume : bool, optional + Run in stochastic volume calculation mode. Defaults to False. + geometry_debug : bool, optional + Turn on geometry debugging during simulation. Defaults to False. + particles : int, optional + Number of particles to simulate per generation. + plot : bool, optional + Run in plotting mode. Defaults to False. + restart_file : str, optional + Path to restart file to use + threads : int, optional + Number of OpenMP threads. If OpenMC is compiled with OpenMP threading + enabled, the default is implementation-dependent but is usually equal + to the number of hardware threads available (or a value set by the + :envvar:`OMP_NUM_THREADS` environment variable). + tracks : bool, optional + Write tracks for all particles. Defaults to False. + event_based : None or bool, optional + Turns on event-based parallelism if True. If None, the value in + the Settings will be used. + openmc_exec : str, optional + Path to OpenMC executable. Defaults to 'openmc'. + mpi_args : list of str, optional + MPI execute command and any additional MPI arguments to pass, + e.g. ['mpiexec', '-n', '8']. + + .. versionadded:: 0.13.0 + + Returns + ------- + args : Iterable of str + The runtime flags converted to CLI arguments of the OpenMC executable + + """ + + args = [openmc_exec] + + if volume: + args.append('--volume') + + if isinstance(particles, Integral) and particles > 0: + args += ['-n', str(particles)] + + if isinstance(threads, Integral) and threads > 0: + args += ['-s', str(threads)] + + if geometry_debug: + args.append('-g') + + if event_based is not None: + if event_based: + args.append('-e') + + if isinstance(restart_file, str): + args += ['-r', restart_file] + + if tracks: + args.append('-t') + + if plot: + args.append('-p') + + if mpi_args is not None: + args = mpi_args + args + + return args def _run(args, output, cwd): @@ -59,12 +137,13 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'): _run([openmc_exec, '-p'], output, cwd) -def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'): +def plot_inline(plots, openmc_exec='openmc', cwd='.'): """Display plots inline in a Jupyter notebook. - This function requires that you have a program installed to convert PPM - files to PNG files. Typically, that would be `ImageMagick - `_ which includes a `convert` command. + .. versionchanged:: 0.13.0 + The *convert_exec* argument was removed since OpenMC now produces + .png images directly. + Parameters ---------- @@ -74,8 +153,6 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'): Path to OpenMC executable cwd : str, optional Path to working directory to run in - convert_exec : str, optional - Command that can convert PPM files into PNG files Raises ------ @@ -83,7 +160,7 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'): If the `openmc` executable returns a non-zero status """ - from IPython.display import Image, display + from IPython.display import display if not isinstance(plots, Iterable): plots = [plots] @@ -94,16 +171,8 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.', convert_exec='convert'): # Run OpenMC in geometry plotting mode plot_geometry(False, openmc_exec, cwd) - images = [] if plots is not None: - for p in plots: - if p.filename is not None: - ppm_file = f'{p.filename}.ppm' - else: - ppm_file = f'plot_{p.id}.ppm' - png_file = ppm_file.replace('.ppm', '.png') - subprocess.check_call([convert_exec, ppm_file, png_file]) - images.append(Image(png_file)) + images = [_get_plot_image(p) for p in plots] display(*images) @@ -126,8 +195,8 @@ def calculate_volumes(threads=None, output=True, cwd='.', ---------- threads : int, optional Number of OpenMP threads. If OpenMC is compiled with OpenMP threading - enabled, the default is implementation-dependent but is usually equal to - the number of hardware threads available (or a value set by the + enabled, the default is implementation-dependent but is usually equal + to the number of hardware threads available (or a value set by the :envvar:`OMP_NUM_THREADS` environment variable). output : bool, optional Capture OpenMC output from standard out @@ -150,12 +219,9 @@ def calculate_volumes(threads=None, output=True, cwd='.', openmc.VolumeCalculation """ - args = [openmc_exec, '--volume'] - if isinstance(threads, Integral) and threads > 0: - args += ['-s', str(threads)] - if mpi_args is not None: - args = mpi_args + args + args = _process_CLI_arguments(volume=True, threads=threads, + openmc_exec=openmc_exec, mpi_args=mpi_args) _run(args, output, cwd) @@ -171,8 +237,8 @@ def run(particles=None, threads=None, geometry_debug=False, Number of particles to simulate per generation. threads : int, optional Number of OpenMP threads. If OpenMC is compiled with OpenMP threading - enabled, the default is implementation-dependent but is usually equal to - the number of hardware threads available (or a value set by the + enabled, the default is implementation-dependent but is usually equal + to the number of hardware threads available (or a value set by the :envvar:`OMP_NUM_THREADS` environment variable). geometry_debug : bool, optional Turn on geometry debugging during simulation. Defaults to False. @@ -201,27 +267,10 @@ def run(particles=None, threads=None, geometry_debug=False, If the `openmc` executable returns a non-zero status """ - args = [openmc_exec] - if isinstance(particles, Integral) and particles > 0: - args += ['-n', str(particles)] - - if isinstance(threads, Integral) and threads > 0: - args += ['-s', str(threads)] - - if geometry_debug: - args.append('-g') - - if event_based: - args.append('-e') - - if isinstance(restart_file, str): - args += ['-r', restart_file] - - if tracks: - args.append('-t') - - if mpi_args is not None: - args = mpi_args + args + args = _process_CLI_arguments( + volume=False, geometry_debug=geometry_debug, particles=particles, + restart_file=restart_file, threads=threads, tracks=tracks, + event_based=event_based, openmc_exec=openmc_exec, mpi_args=mpi_args) _run(args, output, cwd) diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index 82dc92ba4..c14b0d9c2 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -59,3 +59,8 @@ from .tally import * from .settings import settings from .math import * from .plot import * + +# Flag to denote whether or not openmc.lib.init has been called +# TODO: Establish and use a flag in the C++ code to represent the status of the +# openmc_init and openmc_finalize methods +is_initialized = False diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 0da24aab8..e277d527a 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -2,6 +2,7 @@ from contextlib import contextmanager from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p, c_char, POINTER, Structure, c_void_p, create_string_buffer) import sys +import os import numpy as np from numpy.ctypeslib import as_array @@ -109,9 +110,22 @@ def global_bounding_box(): return llc, urc -def calculate_volumes(): - """Run stochastic volume calculation""" - _dll.openmc_calculate_volumes() + +def calculate_volumes(output=True): + """Run stochastic volume calculation + + .. versionchanged:: 0.13.0 + The *output* argument was added. + + Parameters + ---------- + output : bool, optional + Whether or not to show output. Defaults to showing output + + """ + + with quiet_dll(output): + _dll.openmc_calculate_volumes() def current_batch(): @@ -126,13 +140,18 @@ def current_batch(): return c_int.in_dll(_dll, 'current_batch').value -def export_properties(filename=None): +def export_properties(filename=None, output=True): """Export physical properties. + .. versionchanged:: 0.13.0 + The *output* argument was added. + Parameters ---------- filename : str or None Filename to export properties to (defaults to "properties.h5") + output : bool, optional + Whether or not to show output. Defaults to showing output See Also -------- @@ -141,12 +160,15 @@ def export_properties(filename=None): """ if filename is not None: filename = c_char_p(filename.encode()) - _dll.openmc_properties_export(filename) + + with quiet_dll(output): + _dll.openmc_properties_export(filename) def finalize(): """Finalize simulation and free memory""" _dll.openmc_finalize() + openmc.lib.is_initialized = False def find_cell(xyz): @@ -218,15 +240,20 @@ def import_properties(filename): _dll.openmc_properties_import(filename.encode()) -def init(args=None, intracomm=None): +def init(args=None, intracomm=None, output=True): """Initialize OpenMC + .. versionchanged:: 0.13.0 + The *output* argument was added. + Parameters ---------- - args : list of str + args : list of str, optional Command-line arguments - intracomm : mpi4py.MPI.Intracomm or None + intracomm : mpi4py.MPI.Intracomm or None, optional MPI intracommunicator + output : bool, optional + Whether or not to show output. Defaults to showing output """ if args is not None: @@ -252,7 +279,9 @@ def init(args=None, intracomm=None): address = MPI._addressof(intracomm) intracomm = c_void_p(address) - _dll.openmc_init(argc, argv, intracomm) + with quiet_dll(output): + _dll.openmc_init(argc, argv, intracomm) + openmc.lib.is_initialized = True def is_statepoint_batch(): @@ -342,9 +371,20 @@ def next_batch(): return status.value -def plot_geometry(): - """Plot geometry""" - _dll.openmc_plot_geometry() +def plot_geometry(output=True): + """Plot geometry + + .. versionchanged:: 0.13.0 + The *output* argument was added. + + Parameters + ---------- + output : bool, optional + Whether or not to show output. Defaults to showing output + """ + + with quiet_dll(output): + _dll.openmc_plot_geometry() def reset(): @@ -357,9 +397,20 @@ def reset_timers(): _dll.openmc_reset_timers() -def run(): - """Run simulation""" - _dll.openmc_run() +def run(output=True): + """Run simulation + + .. versionchanged:: 0.13.0 + The *output* argument was added. + + Parameters + ---------- + output : bool, optional + Whether or not to show output. Defaults to showing output + """ + + with quiet_dll(output): + _dll.openmc_run() def simulation_init(): @@ -475,3 +526,46 @@ class _FortranObjectWithID(_FortranObject): # assigned. If the array index of the object is out of bounds, an # OutOfBoundsError will be raised here by virtue of referencing self.id self.id + + +@contextmanager +def quiet_dll(output=True): + """This context manager allows us to suppress standard output from DLLs + + Parameters + ---------- + output : bool + Denotes whether the output should be displayed (True) or not (False) + + .. versionadded:: 0.13.0 + + """ + + # This contextmanager is modified from that provided here: + # https://stackoverflow.com/a/14797594 + + if output: + yield + else: + sys.stdout.flush() + # Save the initial file descriptor states + initial_stdout = sys.stdout + initial_stdout_fno = os.dup(sys.stdout.fileno()) + # Get a garbage descriptor so we can throw away output + devnull = os.open(os.devnull, os.O_WRONLY) + + # Get the current stdout stream and make a duplicate of it + new_stdout = os.dup(1) + # Copy the garbage output to the stdout stream + os.dup2(devnull, 1) + os.close(devnull) + # Now point stdout to the re-defined stdout + sys.stdout = os.fdopen(new_stdout, 'w') + + try: + yield + finally: + # Now we just clean up after ourselves and reset the streams + sys.stdout = initial_stdout + sys.stdout.flush() + os.dup2(initial_stdout_fno, 1) diff --git a/openmc/lib/material.py b/openmc/lib/material.py index 6770721ab..fde197d3d 100644 --- a/openmc/lib/material.py +++ b/openmc/lib/material.py @@ -172,7 +172,6 @@ class Material(_FortranObjectWithID): @property def nuclides(self): return self._get_densities()[0] - return nuclides @property def densities(self): diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index 1eab4aa8a..2e9dd18df 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -34,6 +34,7 @@ class _Settings: restart_run = _DLLGlobal(c_bool, 'restart_run') run_CE = _DLLGlobal(c_bool, 'run_CE') verbosity = _DLLGlobal(c_int, 'verbosity') + event_based = _DLLGlobal(c_bool, 'event_based') @property def run_mode(self): diff --git a/openmc/mesh.py b/openmc/mesh.py index 0e35aa6fe..4e1e9565a 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -615,8 +615,8 @@ class UnstructuredMesh(MeshBase): Unique identifier for the mesh name : str Name of the mesh - size : int - Number of elements in the unstructured mesh + length_multiplier: float + Constant multiplier to apply to mesh coordinates Attributes ---------- @@ -626,11 +626,16 @@ class UnstructuredMesh(MeshBase): Name of the mesh filename : str Name of the file containing the unstructured mesh + length_multiplier: float + Multiplicative factor to apply to mesh coordinates library : str Mesh library used for the unstructured mesh tally output : bool Indicates whether or not automatic tally output should be generated for this mesh + specified_length_multiplier: bool + Indicates whether a non-unity length multiplier has been + applied to this mesh volumes : Iterable of float Volumes of the unstructured mesh elements total_volume : float @@ -639,13 +644,16 @@ class UnstructuredMesh(MeshBase): An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0), (1.0, 1.0, 1.0), ...] """ - def __init__(self, filename, library, mesh_id=None, name=''): + 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.library = library self._output = True + self._specified_length_multiplier = False + self.length_multiplier = length_multiplier @property def filename(self): @@ -713,6 +721,20 @@ class UnstructuredMesh(MeshBase): Iterable, Real) self._centroids = centroids + @property + def length_multiplier(self): + return self._length_multiplier + + @length_multiplier.setter + def length_multiplier(self, length_multiplier): + cv.check_type("Unstructured mesh length multiplier", + length_multiplier, + Real) + self._length_multiplier = length_multiplier + + if (self._length_multiplier != 1.0): + self._specified_length_multiplier = True + def __repr__(self): string = super().__repr__() string += '{: <16}=\t{}\n'.format('\tFilename', self.filename) @@ -763,7 +785,7 @@ class UnstructuredMesh(MeshBase): for centroid in self.centroids: # create a point for each centroid - point_id = points.InsertNextPoint(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,)) @@ -817,6 +839,9 @@ class UnstructuredMesh(MeshBase): mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3)) mesh.size = mesh.volumes.size + if 'length_multiplier' in group: + mesh.length_multiplier = group['length_multiplier'][()] + return mesh def to_xml_element(self): @@ -836,6 +861,9 @@ class UnstructuredMesh(MeshBase): subelement = ET.SubElement(element, "filename") subelement.text = self.filename + if (self._specified_length_multiplier): + element.set("length_multiplier", str(self.length_multiplier)) + return element @classmethod @@ -855,5 +883,6 @@ class UnstructuredMesh(MeshBase): mesh_id = int(get_text(elem, 'id')) filename = get_text(elem, 'filename') library = get_text(elem, 'library') + length_multiplier = float(get_text(elem, 'length_multiplier', 1.0)) - return cls(filename, library, mesh_id) + return cls(filename, library, mesh_id, '', length_multiplier) diff --git a/openmc/model/model.py b/openmc/model/model.py index d6c91e73e..bc0796ea7 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,11 +1,32 @@ from collections.abc import Iterable +import operator +import os from pathlib import Path +from numbers import Integral import time +import warnings +import subprocess +from contextlib import contextmanager import h5py import openmc +from openmc.dummy_comm import DummyCommunicator +from openmc.executor import _process_CLI_arguments from openmc.checkvalue import check_type, check_value +from openmc.exceptions import InvalidIDError + + +@contextmanager +def _change_directory(working_dir): + """A context manager for executing in a provided working directory""" + start_dir = Path.cwd() + Path.mkdir(working_dir, exist_ok=True) + os.chdir(working_dir) + try: + yield + finally: + os.chdir(start_dir) class Model: @@ -19,6 +40,10 @@ class Model: not set, it will attempt to create a ``materials.xml`` file based on all materials appearing in the geometry. + .. versionchanged:: 0.13.0 + The model information can now be loaded in to OpenMC directly via + openmc.lib + Parameters ---------- geometry : openmc.Geometry, optional @@ -66,6 +91,30 @@ class Model: if plots is not None: self.plots = plots + # Store dictionaries to the materials and cells by ID and names + if materials is None: + mats = self.geometry.get_all_materials().values() + else: + mats = self.materials + cells = self.geometry.get_all_cells() + # Get the ID maps + self._materials_by_id = {mat.id: mat for mat in mats} + self._cells_by_id = {cell.id: cell for cell in cells.values()} + + # Get the names maps, but since names are not unique, store a list for + # each name key. In this way when the user requests a change by a name, + # the change will be applied to all of the same name. + self._cells_by_name = {} + for cell in cells.values(): + if cell.name not in self._cells_by_name: + self._cells_by_name[cell.name] = set() + self._cells_by_name[cell.name].add(cell) + self._materials_by_name = {} + for mat in mats: + if mat.name not in self._materials_by_name: + self._materials_by_name[mat.name] = set() + self._materials_by_name[mat.name].add(mat) + @property def geometry(self): return self._geometry @@ -86,6 +135,10 @@ class Model: def plots(self): return self._plots + @property + def is_initialized(self): + return openmc.lib.is_initialized + @geometry.setter def geometry(self, geometry): check_type('geometry', geometry, openmc.Geometry) @@ -130,6 +183,8 @@ class Model: def from_xml(cls, geometry='geometry.xml', materials='materials.xml', settings='settings.xml'): """Create model from existing XML files + When initializing this way, the user must manually load plots and + tallies. Parameters ---------- @@ -151,44 +206,159 @@ class Model: settings = openmc.Settings.from_xml(settings) return cls(geometry, materials, settings) - def deplete(self, timesteps, chain_file=None, method='cecm', - fission_q=None, **kwargs): + def init_lib(self, threads=None, geometry_debug=False, restart_file=None, + tracks=False, output=True, event_based=None, intracomm=None): + """Initializes the model in memory via the C API + + .. versionadded:: 0.13.0 + + Parameters + ---------- + threads : int, optional + Number of OpenMP threads. If OpenMC is compiled with OpenMP + threading enabled, the default is implementation-dependent but is + usually equal to the number of hardware threads available + (or a value set by the :envvar:`OMP_NUM_THREADS` environment + variable). + geometry_debug : bool, optional + Turn on geometry debugging during simulation. Defaults to False. + restart_file : str, optional + Path to restart file to use + tracks : bool, optional + Write tracks for all particles. Defaults to False. + output : bool + Capture OpenMC output from standard out + event_based : None or bool, optional + Turns on event-based parallelism if True. If None, the value in + the Settings will be used. + intracomm : mpi4py.MPI.Intracomm or None, optional + MPI intracommunicator + """ + + # TODO: right now the only way to set most of the above parameters via + # the C API are at initialization time despite use-cases existing to + # set them for individual runs. For now this functionality is exposed + # where it exists (here in init), but in the future the functionality + # should be exposed so that it can be accessed via model.run(...) + + args = _process_CLI_arguments( + volume=False, geometry_debug=geometry_debug, + restart_file=restart_file, threads=threads, tracks=tracks, + event_based=event_based) + # Args adds the openmc_exec command in the first entry; remove it + args = args[1:] + + self.finalize_lib() + + # The Model object needs to be aware of the communicator so it can + # use it in certain cases, therefore lets store the communicator + if intracomm is not None: + self._intracomm = intracomm + else: + self._intracomm = DummyCommunicator() + + if self._intracomm.rank == 0: + self.export_to_xml() + self._intracomm.barrier() + + # We cannot pass DummyCommunicator to openmc.lib.init so pass instead + # the user-provided intracomm which will either be None or an mpi4py + # communicator + openmc.lib.init(args=args, intracomm=intracomm, output=output) + + def finalize_lib(self): + """Finalize simulation and free memory allocated for the C API + + .. versionadded:: 0.13.0 + + """ + + openmc.lib.finalize() + + def deplete(self, timesteps, method='cecm', final_step=True, + operator_kwargs=None, directory='.', output=True, + **integrator_kwargs): """Deplete model using specified timesteps/power + .. versionchanged:: 0.13.0 + The *final_step*, *operator_kwargs*, *directory*, and *output* + arguments were added. + Parameters ---------- timesteps : iterable of float Array of timesteps in units of [s]. Note that values are not cumulative. - chain_file : str, optional - Path to the depletion chain XML file. Defaults to the chain - found under the ``depletion_chain`` in the - :envvar:`OPENMC_CROSS_SECTIONS` environment variable if it exists. - method : str - Integration method used for depletion (e.g., 'cecm', 'predictor') - fission_q : dict, optional - Dictionary of nuclides and their fission Q values [eV]. - If not given, values will be pulled from the ``chain_file``. - **kwargs - Keyword arguments passed to integration function (e.g., - :func:`openmc.deplete.integrator.cecm`) + method : str, optional + Integration method used for depletion (e.g., 'cecm', 'predictor'). + Defaults to 'cecm'. + final_step : bool, optional + Indicate whether or not a transport solve should be run at the end + of the last timestep. Defaults to running this transport solve. + operator_kwargs : dict + Keyword arguments passed to the depletion Operator initializer + (e.g., :func:`openmc.deplete.Operator`) + directory : str, optional + Directory to write XML files to. If it doesn't exist already, it + will be created. Defaults to the current working directory + output : bool + Capture OpenMC output from standard out + integrator_kwargs : dict + Remaining keyword arguments passed to the depletion Integrator + initializer (e.g., :func:`openmc.deplete.integrator.cecm`). """ - # Import the depletion module. This is done here rather than the module - # header to delay importing openmc.lib (through openmc.deplete) which - # can be tough to install properly. + + if operator_kwargs is None: + op_kwargs = {} + elif isinstance(operator_kwargs, dict): + op_kwargs = operator_kwargs + else: + raise ValueError("operator_kwargs must be a dict or None") + + # Import openmc.deplete here so the Model can be used even if the + # shared library is unavailable. import openmc.deplete as dep - # Create OpenMC transport operator - op = dep.Operator( - self.geometry, self.settings, chain_file, - fission_q=fission_q, - ) + # Store whether or not the library was initialized when we started + started_initialized = self.is_initialized - # Perform depletion - check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4', - 'si_celi', 'si_leqi', 'celi', 'leqi')) - getattr(dep.integrator, method)(op, timesteps, **kwargs) + with _change_directory(Path(directory)): + with openmc.lib.quiet_dll(output): + depletion_operator = \ + dep.Operator(self.geometry, self.settings, **op_kwargs) + + # Tell depletion_operator.finalize NOT to clear C API memory when + # it is done + depletion_operator.cleanup_when_done = False + + # Set up the integrator + check_value('method', method, + dep.integrators.integrator_by_name.keys()) + integrator_class = dep.integrators.integrator_by_name[method] + integrator = integrator_class(depletion_operator, timesteps, + **integrator_kwargs) + + # Now perform the depletion + with openmc.lib.quiet_dll(output): + integrator.integrate(final_step) + + # Now make the python Materials match the C API material data + for mat_id, mat in self._materials_by_id.items(): + if mat.depletable: + # Get the C data + c_mat = openmc.lib.materials[mat_id] + nuclides, densities = c_mat._get_densities() + # And now we can remove isotopes and add these ones in + mat.nuclides.clear() + for nuc, density in zip(nuclides, densities): + mat.add_nuclide(nuc, density) + mat.set_density('atom/b-cm', sum(densities)) + + # If we didnt start intialized, we should cleanup after ourselves + if not started_initialized: + depletion_operator.cleanup_when_done = True + depletion_operator.finalize() def export_to_xml(self, directory='.'): """Export model to XML files. @@ -226,6 +396,9 @@ class Model: def import_properties(self, filename): """Import physical properties + .. versionchanged:: 0.13.0 + This method now updates values as loaded in memory with the C API + Parameters ---------- filename : str @@ -254,32 +427,69 @@ class Model: cell = cells[cell_id] if cell.fill_type in ('material', 'distribmat'): cell.temperature = group['temperature'][()] + if self.is_initialized: + lib_cell = openmc.lib.cells[cell_id] + lib_cell.set_temperature(group['temperature'][()]) # Make sure number of materials matches mats_group = fh['materials'] n_cells = mats_group.attrs['n_materials'] if n_cells != len(materials): - raise ValueError("Number of materials in properties file doesn't " - "match current model.") + raise ValueError("Number of materials in properties file " + "doesn't match current model.") # Update material densities for name, group in mats_group.items(): mat_id = int(name.split()[1]) atom_density = group.attrs['atom_density'] materials[mat_id].set_density('atom/b-cm', atom_density) + if self.is_initialized: + C_mat = openmc.lib.materials[mat_id] + C_mat.set_density(atom_density, 'atom/b-cm') - def run(self, **kwargs): - """Creates the XML files, runs OpenMC, and returns the path to the last + def run(self, particles=None, threads=None, geometry_debug=False, + restart_file=None, tracks=False, output=True, cwd='.', + openmc_exec='openmc', mpi_args=None, event_based=None): + """Runs OpenMC. If the C API has been initialized, then the C API is + used, otherwise, this method creates the XML files and runs OpenMC via + a system call. In both cases this method returns the path to the last statepoint file generated. .. versionchanged:: 0.12 Instead of returning the final k-effective value, this function now returns the path to the final statepoint written. + .. versionchanged:: 0.13.0 + This method can utilize the C API for execution + Parameters ---------- - **kwargs - Keyword arguments passed to :func:`openmc.run` + particles : int, optional + Number of particles to simulate per generation. + threads : int, optional + Number of OpenMP threads. If OpenMC is compiled with OpenMP + threading enabled, the default is implementation-dependent but is + usually equal to the number of hardware threads available (or a + value set by the :envvar:`OMP_NUM_THREADS` environment variable). + geometry_debug : bool, optional + Turn on geometry debugging during simulation. Defaults to False. + restart_file : str, optional + Path to restart file to use + tracks : bool, optional + Write tracks for all particles. Defaults to False. + output : bool, optional + Capture OpenMC output from standard out + cwd : str, optional + Path to working directory to run in. Defaults to the current + working directory. + openmc_exec : str, optional + Path to OpenMC executable. Defaults to 'openmc'. + mpi_args : list of str, optional + MPI execute command and any additional MPI arguments to pass, + e.g. ['mpiexec', '-n', '8']. + event_based : None or bool, optional + Turns on event-based parallelism if True. If None, the value in + the Settings will be used. Returns ------- @@ -289,23 +499,346 @@ class Model: """ - self.export_to_xml() - - # Setting tstart here ensures we don't pick up any pre-existing statepoint - # files in the output directory + # Setting tstart here ensures we don't pick up any pre-existing + # statepoint files in the output directory tstart = time.time() last_statepoint = None - openmc.run(**kwargs) + # Operate in the provided working directory + with _change_directory(Path(cwd)): + if self.is_initialized: + # Handle the run options as applicable + # First dont allow ones that must be set via init + for arg_name, arg, default in zip( + ['threads', 'geometry_debug', 'restart_file', 'tracks'], + [threads, geometry_debug, restart_file, tracks], + [None, False, None, False]): + if arg != default: + msg = f"{arg_name} must be set via Model.is_initialized(...)" + raise ValueError(msg) - # Get output directory and return the last statepoint written by this run - if self.settings.output and 'path' in self.settings.output: - output_dir = Path(self.settings.output['path']) - else: - output_dir = Path.cwd() - for sp in output_dir.glob('statepoint.*.h5'): - mtime = sp.stat().st_mtime - if mtime >= tstart: # >= allows for poor clock resolution - tstart = mtime - last_statepoint = sp + init_particles = openmc.lib.settings.particles + if particles is not None: + if isinstance(particles, Integral) and particles > 0: + openmc.lib.settings.particles = particles + + init_event_based = openmc.lib.settings.event_based + if event_based is not None: + openmc.lib.settings.event_based = event_based + + # Then run using the C API + openmc.lib.run(output) + + # Reset changes for the openmc.run kwargs handling + openmc.lib.settings.particles = init_particles + openmc.lib.settings.event_based = init_event_based + + else: + # Then run via the command line + self.export_to_xml() + openmc.run(particles, threads, geometry_debug, restart_file, + tracks, output, Path('.'), openmc_exec, mpi_args, + event_based) + + # Get output directory and return the last statepoint written + if self.settings.output and 'path' in self.settings.output: + output_dir = Path(self.settings.output['path']) + else: + output_dir = Path.cwd() + for sp in output_dir.glob('statepoint.*.h5'): + mtime = sp.stat().st_mtime + if mtime >= tstart: # >= allows for poor clock resolution + tstart = mtime + last_statepoint = sp return last_statepoint + + def calculate_volumes(self, threads=None, output=True, cwd='.', + openmc_exec='openmc', mpi_args=None, + apply_volumes=True): + """Runs an OpenMC stochastic volume calculation and, if requested, + applies volumes to the model + + .. versionadded:: 0.13.0 + + Parameters + ---------- + threads : int, optional + Number of OpenMP threads. If OpenMC is compiled with OpenMP + threading enabled, the default is implementation-dependent but is + usually equal to the number of hardware threads available (or a + value set by the :envvar:`OMP_NUM_THREADS` environment variable). + This currenty only applies to the case when not using the C API. + output : bool, optional + Capture OpenMC output from standard out + openmc_exec : str, optional + Path to OpenMC executable. Defaults to 'openmc'. + This only applies to the case when not using the C API. + mpi_args : list of str, optional + MPI execute command and any additional MPI arguments to pass, + e.g. ['mpiexec', '-n', '8']. + This only applies to the case when not using the C API. + cwd : str, optional + Path to working directory to run in. Defaults to the current + working directory. + apply_volumes : bool, optional + Whether apply the volume calculation results from this calculation + to the model. Defaults to applying the volumes. + """ + + if len(self.settings.volume_calculations) == 0: + # Then there is no volume calculation specified + raise ValueError("The Settings.volume_calculation attribute must" + " be specified before executing this method!") + + with _change_directory(Path(cwd)): + if self.is_initialized: + if threads is not None: + msg = "Threads must be set via Model.is_initialized(...)" + raise ValueError(msg) + if mpi_args is not None: + msg = "The MPI environment must be set otherwise such as" \ + "with the call to mpi4py" + raise ValueError(msg) + + # Compute the volumes + openmc.lib.calculate_volumes(output) + + else: + self.export_to_xml() + openmc.calculate_volumes(threads=threads, output=output, + openmc_exec=openmc_exec, + mpi_args=mpi_args) + + # Now we apply the volumes + if apply_volumes: + # Load the results and add them to the model + for i, vol_calc in enumerate(self.settings.volume_calculations): + vol_calc.load_results(f"volume_{i + 1}.h5") + # First add them to the Python side + self.geometry.add_volume_information(vol_calc) + + # And now repeat for the C API + if self.is_initialized and vol_calc.domain_type == 'material': + # Then we can do this in the C API + for domain_id in vol_calc.ids: + openmc.lib.materials[domain_id].volume = \ + vol_calc.volumes[domain_id].n + + def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'): + """Creates plot images as specified by the Model.plots attribute + + .. versionadded:: 0.13.0 + + Parameters + ---------- + output : bool, optional + Capture OpenMC output from standard out + cwd : str, optional + Path to working directory to run in. Defaults to the current + working directory. + openmc_exec : str, optional + Path to OpenMC executable. Defaults to 'openmc'. + This only applies to the case when not using the C API. + + """ + + if len(self.plots) == 0: + # Then there is no volume calculation specified + raise ValueError("The Model.plots attribute must be specified " + "before executing this method!") + + with _change_directory(Path(cwd)): + if self.is_initialized: + # Compute the volumes + openmc.lib.plot_geometry(output) + else: + self.export_to_xml() + openmc.plot_geometry(output=output, openmc_exec=openmc_exec) + + def _change_py_lib_attribs(self, names_or_ids, value, obj_type, + attrib_name, density_units='atom/b-cm'): + # Method to do the same work whether it is a cell or material and + # a temperature or volume + check_type('names_or_ids', names_or_ids, Iterable, (Integral, str)) + check_type('obj_type', obj_type, str) + obj_type = obj_type.lower() + check_value('obj_type', obj_type, ('material', 'cell')) + check_value('attrib_name', attrib_name, + ('temperature', 'volume', 'density', 'rotation', + 'translation')) + # The C API only allows setting density units of atom/b-cm and g/cm3 + check_value('density_units', density_units, ('atom/b-cm', 'g/cm3')) + # The C API has no way to set cell volume or material temperature + # so lets raise exceptions as needed + if obj_type == 'cell' and attrib_name == 'volume': + raise NotImplementedError( + 'Setting a Cell volume is not supported!') + if obj_type == 'material' and attrib_name == 'temperature': + raise NotImplementedError( + 'Setting a material temperature is not supported!') + + # And some items just dont make sense + if obj_type == 'cell' and attrib_name == 'density': + raise ValueError('Cannot set a Cell density!') + if obj_type == 'material' and attrib_name in ('rotation', + 'translation'): + raise ValueError('Cannot set a material rotation/translation!') + + # Set the + if obj_type == 'cell': + by_name = self._cells_by_name + by_id = self._cells_by_id + obj_by_id = openmc.lib.cells + else: + by_name = self._materials_by_name + by_id = self._materials_by_id + obj_by_id = openmc.lib.materials + + # Get the list of ids to use if converting from names and accepting + # only values that have actual ids + ids = [] + for name_or_id in names_or_ids: + if isinstance(name_or_id, Integral): + if name_or_id in by_id: + ids.append(int(name_or_id)) + else: + cap_obj = obj_type.capitalize() + msg = f'{cap_obj} ID {name_or_id} " \ + "is not present in the model!' + raise InvalidIDError(msg) + elif isinstance(name_or_id, str): + if name_or_id in by_name: + # Then by_name[name_or_id] is a list so we need to add all + # entries + ids.extend([obj.id for obj in by_name[name_or_id]]) + else: + cap_obj = obj_type.capitalize() + msg = f'{cap_obj} {name_or_id} " \ + "is not present in the model!' + raise InvalidIDError(msg) + + # Now perform the change to both python and C API + for id_ in ids: + obj = by_id[id_] + if attrib_name == 'density': + obj.set_density(density_units, value) + else: + setattr(obj, attrib_name, value) + # Next lets keep what is in C API memory up to date as well + if self.is_initialized: + lib_obj = obj_by_id[id_] + if attrib_name == 'density': + lib_obj.set_density(value, density_units) + elif attrib_name == 'temperature': + lib_obj.set_temperature(value) + else: + setattr(lib_obj, attrib_name, value) + + def rotate_cells(self, names_or_ids, vector): + """Rotate the identified cell(s) by the specified rotation vector. + The rotation is only applied to cells filled with a universe. + + .. note:: If applying this change to a name that is not unique, then + the change will be applied to all objects of that name. + + .. versionadded:: 0.13.0 + + Parameters + ---------- + names_or_ids : Iterable of str or int + The cell names (if str) or id (if int) that are to be translated + or rotated. This parameter can include a mix of names and ids. + vector : Iterable of float + The rotation vector of length 3 to apply. This array specifies the + angles in degrees about the x, y, and z axes, respectively. + + """ + + self._change_py_lib_attribs(names_or_ids, vector, 'cell', 'rotation') + + def translate_cells(self, names_or_ids, vector): + """Translate the identified cell(s) by the specified translation vector. + The translation is only applied to cells filled with a universe. + + .. note:: If applying this change to a name that is not unique, then + the change will be applied to all objects of that name. + + .. versionadded:: 0.13.0 + + Parameters + ---------- + names_or_ids : Iterable of str or int + The cell names (if str) or id (if int) that are to be translated + or rotated. This parameter can include a mix of names and ids. + vector : Iterable of float + The translation vector of length 3 to apply. This array specifies + the x, y, and z dimensions of the translation. + + """ + + self._change_py_lib_attribs(names_or_ids, vector, 'cell', + 'translation') + + def update_densities(self, names_or_ids, density, density_units='atom/b-cm'): + """Update the density of a given set of materials to a new value + + .. note:: If applying this change to a name that is not unique, then + the change will be applied to all objects of that name. + + .. versionadded:: 0.13.0 + + Parameters + ---------- + names_or_ids : Iterable of str or int + The material names (if str) or id (if int) that are to be updated. + This parameter can include a mix of names and ids. + density : float + The density to apply in the units specified by `density_units` + density_units : {'atom/b-cm', 'g/cm3'}, optional + Units for `density`. Defaults to 'atom/b-cm' + + """ + + self._change_py_lib_attribs(names_or_ids, density, 'material', + 'density', density_units) + + def update_cell_temperatures(self, names_or_ids, temperature): + """Update the temperature of a set of cells to the given value + + .. note:: If applying this change to a name that is not unique, then + the change will be applied to all objects of that name. + + .. versionadded:: 0.13.0 + + Parameters + ---------- + names_or_ids : Iterable of str or int + The cell names (if str) or id (if int) that are to be updated. + This parameter can include a mix of names and ids. + temperature : float + The temperature to apply in units of Kelvin + + """ + + self._change_py_lib_attribs(names_or_ids, temperature, 'cell', + 'temperature') + + def update_material_volumes(self, names_or_ids, volume): + """Update the volume of a set of materials to the given value + + .. note:: If applying this change to a name that is not unique, then + the change will be applied to all objects of that name. + + .. versionadded:: 0.13.0 + + Parameters + ---------- + names_or_ids : Iterable of str or int + The material names (if str) or id (if int) that are to be updated. + This parameter can include a mix of names and ids. + volume : float + The volume to apply in units of cm^3 + + """ + + self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume') diff --git a/openmc/mpi.py b/openmc/mpi.py new file mode 100644 index 000000000..cfc10c0d1 --- /dev/null +++ b/openmc/mpi.py @@ -0,0 +1,8 @@ +try: + from mpi4py import MPI + comm = MPI.COMM_WORLD +except ImportError: + from unittest.mock import Mock + MPI = Mock() + from openmc.dummy_comm import DummyCommunicator + comm = DummyCommunicator() diff --git a/openmc/plots.py b/openmc/plots.py index 0abff20d7..521be95ad 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -1,6 +1,7 @@ from collections.abc import Iterable, Mapping from numbers import Real, Integral from pathlib import Path +import shutil import subprocess from xml.etree import ElementTree as ET @@ -165,6 +166,18 @@ _SVG_COLORS = { } +def _get_plot_image(plot): + from IPython.display import Image + + # Make sure .png file was created + stem = plot.filename if plot.filename is not None else f'plot_{plot.id}' + png_file = f'{stem}.png' + if not Path(png_file).exists(): + raise FileNotFoundError(f"Could not find .png image for plot {plot.id}") + + return Image(png_file) + + class Plot(IDManagerMixin): """Definition of a finite region of space to be plotted. @@ -671,15 +684,14 @@ class Plot(IDManagerMixin): return element - def to_ipython_image(self, openmc_exec='openmc', cwd='.', - convert_exec='convert'): + def to_ipython_image(self, openmc_exec='openmc', cwd='.'): """Render plot as an image - This method runs OpenMC in plotting mode to produce a bitmap image which - is then converted to a .png file and loaded in as an - :class:`IPython.display.Image` object. As such, it requires that your - model geometry, materials, and settings have already been exported to - XML. + This method runs OpenMC in plotting mode to produce a .png file. + + .. versionchanged:: 0.13.0 + The *convert_exec* argument was removed since OpenMC now produces + .png images directly. Parameters ---------- @@ -687,8 +699,6 @@ class Plot(IDManagerMixin): Path to OpenMC executable cwd : str, optional Path to working directory to run in - convert_exec : str, optional - Command that can convert PPM files into PNG files Returns ------- @@ -696,23 +706,14 @@ class Plot(IDManagerMixin): Image generated """ - from IPython.display import Image - # Create plots.xml Plots([self]).export_to_xml() # Run OpenMC in geometry plotting mode openmc.plot_geometry(False, openmc_exec, cwd) - # Convert to .png - if self.filename is not None: - ppm_file = f'{self.filename}.ppm' - else: - ppm_file = f'plot_{self.id}.ppm' - png_file = ppm_file.replace('.ppm', '.png') - subprocess.check_call([convert_exec, ppm_file, png_file]) - - return Image(png_file) + # Return produced image + return _get_plot_image(self) class Plots(cv.CheckedList): diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 4bbc50927..af6c9a6d1 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -812,8 +812,48 @@ class Mixture(Univariate): self._distribution = distribution def to_xml_element(self, element_name): - raise NotImplementedError + """Return XML representation of the mixture distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing mixture distribution data + + """ + element = ET.Element(element_name) + element.set("type", "mixture") + + for p, d in zip(self.probability, self.distribution): + data = ET.SubElement(element, "pair") + data.set("probability", str(p)) + data.append(d.to_xml_element("dist")) + + return element @classmethod def from_xml_element(cls, elem): - raise NotImplementedError + """Generate mixture distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Mixture + Mixture distribution generated from XML element + + """ + probability = [] + distribution = [] + for pair in elem.findall('pair'): + probability.append(float(get_text(pair, 'probability'))) + distribution.append(Univariate.from_xml_element(pair.find("dist"))) + + return cls(probability, distribution) diff --git a/src/cell.cpp b/src/cell.cpp index 109b66dde..fe138a089 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -10,7 +10,7 @@ #include #include -#include +#include #include "openmc/capi.h" #include "openmc/constants.h" diff --git a/src/dagmc.cpp b/src/dagmc.cpp index a6ff0e3b8..9b894d3bb 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -713,7 +713,13 @@ int32_t next_cell( namespace openmc { -void read_dagmc_universes(pugi::xml_node node) {}; +void read_dagmc_universes(pugi::xml_node node) +{ + if (check_for_node(node, "dagmc_universe")) { + fatal_error("DAGMC Universes are present but OpenMC was not configured" + "with DAGMC"); + } +}; void check_dagmc_root_univ() {}; } // namespace openmc diff --git a/src/distribution.cpp b/src/distribution.cpp index 4713110ab..62296cb25 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -7,6 +7,8 @@ #include // for runtime_error #include // for string, stod +#include + #include "openmc/error.h" #include "openmc/math_functions.h" #include "openmc/random_dist.h" @@ -287,6 +289,48 @@ double Equiprobable::sample(uint64_t* seed) const return xl + ((n - 1) * r - i) * (xr - xl); } +//============================================================================== +// Mixture implementation +//============================================================================== + +Mixture::Mixture(pugi::xml_node node) +{ + double cumsum = 0.0; + for (pugi::xml_node pair : node.children("pair")) { + // Check that required data exists + if (!pair.attribute("probability")) fatal_error("Mixture pair element does not have probability."); + if (!pair.child("dist")) fatal_error("Mixture pair element does not have a distribution."); + + // cummulative sum of probybilities + cumsum += std::stod(pair.attribute("probability").value()); + + // Save cummulative probybility and distrubution + distribution_.push_back( + std::make_pair(cumsum, distribution_from_xml(pair.child("dist")))); + } + + // Normalize cummulative probabilities to 1 + for (auto& pair : distribution_) { + pair.first /= cumsum; + } +} + +double Mixture::sample(uint64_t* seed) const +{ + // Sample value of CDF + const double p = prn(seed); + + // find matching distribution + const auto it = std::lower_bound(distribution_.cbegin(), distribution_.cend(), + p, [](const DistPair& pair, double p) { return pair.first < p; }); + + // This should not happen. Catch it + Ensures(it != distribution_.cend()); + + // Sample the chosen distribution + return it->second->sample(seed); +} + //============================================================================== // Helper function //============================================================================== @@ -315,6 +359,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node) dist = UPtrDist {new Discrete(node)}; } else if (type == "tabular") { dist = UPtrDist {new Tabular(node)}; + } else if (type == "mixture") { + dist = UPtrDist {new Mixture(node)}; } else { openmc::fatal_error("Invalid distribution type: " + type); } diff --git a/src/finalize.cpp b/src/finalize.cpp index b1e09a37e..2f6c65aaf 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -15,6 +15,7 @@ #include "openmc/message_passing.h" #include "openmc/nuclide.h" #include "openmc/photon.h" +#include "openmc/plot.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" #include "openmc/simulation.h" @@ -45,6 +46,7 @@ void free_memory() free_memory_mesh(); free_memory_tally(); free_memory_bank(); + free_memory_plot(); if (mpi::master) { free_memory_cmfd(); } @@ -128,6 +130,7 @@ int openmc_finalize() data::temperature_min = 0.0; data::temperature_max = INFTY; model::root_universe = -1; + model::plotter_seed = 1; openmc::openmc_set_seed(DEFAULT_SEED); // Deallocate arrays diff --git a/src/initialize.cpp b/src/initialize.cpp index e53767691..75ec20518 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -305,9 +305,11 @@ void read_input_xml() // Initialize distribcell_filters prepare_distribcell(); + // Read the plots.xml regardless of plot mode in case plots are requested + // via the API + read_plots_xml(); if (settings::run_mode == RunMode::PLOTTING) { // Read plots.xml if it exists - read_plots_xml(); if (mpi::master && settings::verbosity >= 5) print_plot(); diff --git a/src/mesh.cpp b/src/mesh.cpp index 4798bfcf5..1dc7d3402 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -2,7 +2,7 @@ #include // for copy, equal, min, min_element #include // for ceil #include // for size_t -#include +#include #include #ifdef OPENMC_MPI @@ -35,6 +35,7 @@ #ifdef LIBMESH #include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" +#include "libmesh/mesh_modification.h" #endif namespace openmc { @@ -171,6 +172,12 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) } } + // check if a length unit multiplier was specified + if (check_for_node(node, "length_multiplier")) { + length_multiplier_ = std::stod(get_node_value(node, "length_multiplier")); + specified_length_multiplier_ = true; + } + // get the filename of the unstructured mesh to load if (check_for_node(node, "filename")) { filename_ = get_node_value(node, "filename"); @@ -218,9 +225,21 @@ void UnstructuredMesh::to_hdf5(hid_t group) const 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_); + close_group(mesh_group); } +void UnstructuredMesh::set_length_multiplier(double length_multiplier) +{ + length_multiplier_ = length_multiplier; + + if (length_multiplier_ != 1.0) + specified_length_multiplier_ = true; +} + void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const { *in_mesh = true; @@ -1585,9 +1604,10 @@ MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node) initialize(); } -MOABMesh::MOABMesh(const std::string& filename) +MOABMesh::MOABMesh(const std::string& filename, double length_multiplier) { filename_ = filename; + set_length_multiplier(length_multiplier); initialize(); } @@ -1634,6 +1654,34 @@ void MOABMesh::initialize() fatal_error("Failed to add tetrahedra to an entity set."); } + if (specified_length_multiplier_) { + // get the connectivity of all tets + moab::Range adj; + rval = mbi_->get_adjacencies(ehs_, 0, true, adj, moab::Interface::UNION); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to get adjacent vertices of tetrahedra."); + } + // scale all vertex coords by multiplier (done individually so not all + // coordinates are in memory twice at once) + for (auto vert : adj) { + // retrieve coords + std::array coord; + rval = mbi_->get_coords(&vert, 1, coord.data()); + if (rval != moab::MB_SUCCESS) { + fatal_error("Could not get coordinates of vertex."); + } + // scale coords + for (auto& c : coord) { + c *= length_multiplier_; + } + // set new coords + rval = mbi_->set_coords(&vert, 1, coord.data()); + if (rval != moab::MB_SUCCESS) { + fatal_error("Failed to set new vertex coordinates"); + } + } + } + // build acceleration data structures compute_barycentric_data(ehs_); build_kdtree(ehs_); @@ -2144,9 +2192,10 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) initialize(); } -LibMesh::LibMesh(const std::string& filename) +LibMesh::LibMesh(const std::string& filename, double length_multiplier) { filename_ = filename; + set_length_multiplier(length_multiplier); initialize(); } @@ -2162,6 +2211,11 @@ void LibMesh::initialize() 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(); // ensure that the loaded mesh is 3 dimensional diff --git a/src/output.cpp b/src/output.cpp index 73f3a0e9a..7bc012296 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -72,26 +72,26 @@ void title() // Write version information fmt::print( - " | The OpenMC Monte Carlo Code\n" - " Copyright | 2011-2021 MIT and OpenMC contributors\n" - " License | https://docs.openmc.org/en/latest/license.html\n" - " Version | {}.{}.{}{}\n", + " | The OpenMC Monte Carlo Code\n" + " Copyright | 2011-2021 MIT, UChicago Argonne LLC, and contributors\n" + " License | https://docs.openmc.org/en/latest/license.html\n" + " Version | {}.{}.{}{}\n", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE, VERSION_DEV ? "-dev" : ""); #ifdef GIT_SHA1 - fmt::print(" Git SHA1 | {}\n", GIT_SHA1); + fmt::print(" Git SHA1 | {}\n", GIT_SHA1); #endif // Write the date and time - fmt::print(" Date/Time | {}\n", time_stamp()); + fmt::print(" Date/Time | {}\n", time_stamp()); #ifdef OPENMC_MPI // Write number of processors - fmt::print(" MPI Processes | {}\n", mpi::n_procs); + fmt::print(" MPI Processes | {}\n", mpi::n_procs); #endif #ifdef _OPENMP // Write number of OpenMP threads - fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); + fmt::print(" OpenMP Threads | {}\n", omp_get_max_threads()); #endif std::cout << std::endl; } @@ -335,8 +335,8 @@ void print_version() #ifdef GIT_SHA1 fmt::print("Git SHA1: {}\n", GIT_SHA1); #endif - fmt::print("Copyright (c) 2011-2021 Massachusetts Institute of " - "Technology and OpenMC contributors\nMIT/X license at " + fmt::print("Copyright (c) 2011-2021 MIT, UChicago Argonne LLC, and " + "contributors\nMIT/X license at " "\n"); } } diff --git a/src/plot.cpp b/src/plot.cpp index 789437691..918300e09 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1,12 +1,16 @@ #include "openmc/plot.h" #include +#include #include #include #include "xtensor/xview.hpp" #include #include +#ifdef USE_LIBPNG +#include +#endif #include "openmc/constants.h" #include "openmc/error.h" @@ -103,25 +107,29 @@ uint64_t plotter_seed = 1; extern "C" int openmc_plot_geometry() { + for (auto& pl : model::plots) { write_message(5, "Processing plot {}: {}...", pl.id_, pl.path_plot_); if (PlotType::slice == pl.type_) { // create 2D image - create_ppm(pl); + create_image(pl); } else if (PlotType::voxel == pl.type_) { // create voxel file for 3D viewing create_voxel(pl); } } + return 0; } void read_plots_xml() { - // Check if plots.xml exists + // Check if plots.xml exists; this is only necessary when the plot runmode is + // initiated. Otherwise, we want to read plots.xml because it may be called + // later via the API. In that case, its ok for a plots.xml to not exist std::string filename = settings::path_input + "plots.xml"; - if (!file_exists(filename)) { + if (!file_exists(filename) && settings::run_mode == RunMode::PLOTTING) { fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename)); } @@ -138,12 +146,18 @@ void read_plots_xml() } } +void free_memory_plot() +{ + model::plots.clear(); + model::plot_map.clear(); +} + //============================================================================== -// CREATE_PPM creates an image based on user input from a plots.xml -// specification in the portable pixmap format (PPM) +// CREATE_IMAGE creates an image based on user input from a plots.xml +// specification in the PNG/PPM format //============================================================================== -void create_ppm(Plot const& pl) +void create_image(Plot const& pl) { size_t width = pl.pixels_[0]; @@ -184,8 +198,12 @@ void create_ppm(Plot const& pl) draw_mesh_lines(pl, data); } - // write ppm data to file +// create image file +#ifdef USE_LIBPNG + output_png(pl, data); +#else output_ppm(pl, data); +#endif } void Plot::set_id(pugi::xml_node plot_node) @@ -238,7 +256,11 @@ void Plot::set_output_path(pugi::xml_node plot_node) // add appropriate file extension to name switch (type_) { case PlotType::slice: +#ifdef USE_LIBPNG + filename.append(".png"); +#else filename.append(".ppm"); +#endif break; case PlotType::voxel: filename.append(".h5"); @@ -673,6 +695,60 @@ void output_ppm(Plot const& pl, const ImageData& data) of << "\n"; } +//============================================================================== +// OUTPUT_PNG writes out a previously generated image to a PNG file +//============================================================================== + +#ifdef USE_LIBPNG +void output_png(Plot const& pl, const ImageData& data) +{ + // Open PNG file for writing + std::string fname = pl.path_plot_; + fname = strtrim(fname); + auto fp = std::fopen(fname.c_str(), "wb"); + + // Initialize write and info structures + auto png_ptr = + png_create_write_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr); + auto info_ptr = png_create_info_struct(png_ptr); + + // Setup exception handling + if (setjmp(png_jmpbuf(png_ptr))) + fatal_error("Error during png creation"); + + png_init_io(png_ptr, fp); + + // Write header (8 bit colour depth) + int width = pl.pixels_[0]; + int height = pl.pixels_[1]; + png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB, + PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); + png_write_info(png_ptr, info_ptr); + + // Allocate memory for one row (3 bytes per pixel - RGB) + std::vector row(3 * width); + + // Write color for each pixel + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + RGBColor rgb = data(x, y); + row[3 * x] = rgb.red; + row[3 * x + 1] = rgb.green; + row[3 * x + 2] = rgb.blue; + } + png_write_row(png_ptr, row.data()); + } + + // End write + png_write_end(png_ptr, nullptr); + + // Clean up data structures + std::fclose(fp); + png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1); + png_destroy_write_struct(&png_ptr, nullptr); +} +#endif + //============================================================================== // DRAW_MESH_LINES draws mesh line boundaries on an image //============================================================================== @@ -777,7 +853,7 @@ void draw_mesh_lines(Plot const& pl, ImageData& data) //============================================================================== // CREATE_VOXEL outputs a binary file that can be input into silomesh for 3D -// geometry visualization. It works the same way as create_ppm by dragging a +// geometry visualization. It works the same way as create_image by dragging a // particle across the geometry for the specified number of voxels. The first 3 // int's in the binary are the number of x, y, and z voxels. The next 3 // double's are the widths of the voxels in the x, y, and z directions. The diff --git a/src/surface.cpp b/src/surface.cpp index c19de1e00..73a5aaaf0 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -5,7 +5,7 @@ #include #include -#include +#include #include "openmc/array.h" #include "openmc/container_util.h" diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index a5f52e95a..3f895b4f8 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -1,7 +1,7 @@ #include "openmc/tallies/filter_mesh.h" #include -#include +#include #include "openmc/capi.h" #include "openmc/constants.h" diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index 29976f3a1..359df379b 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -3,7 +3,7 @@ #include // For pair #include -#include +#include #include "openmc/capi.h" #include "openmc/error.h" diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index 65004983c..eb4c8bdfd 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -5,7 +5,7 @@ #include // For pair #include -#include +#include #include "openmc/capi.h" #include "openmc/error.h" diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index 54b6a253f..9174202fb 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -a8192f6029cf99748816fab2618fa48e180656eba313940ccdfff3e56890a5dadd13bf92cd7e3be6a4b535bca5e2a58a905fcfba018fb922273f8be6dfc19859 \ No newline at end of file +f85c20735a0c08525fe48b19a8e075c074539ee6c8860268fa0cb515d842496b7086e5b94305ef78dcf2106bb193abdf438259ac8ff1d0245a2782eb6f5af873 \ No newline at end of file diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index a6f7306a0..2fd6a5b79 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -19,7 +19,7 @@ class PlotTestHarness(TestHarness): openmc.plot_geometry(openmc_exec=config['exe']) def _test_output_created(self): - """Make sure *.ppm has been created.""" + """Make sure *.png has been created.""" for fname in self._plot_names: assert os.path.exists(fname), 'Plot output file does not exist.' @@ -34,8 +34,8 @@ class PlotTestHarness(TestHarness): outstr = bytes() for fname in self._plot_names: - if fname.endswith('.ppm'): - # Add PPM output to results + if fname.endswith('.png'): + # Add PNG output to results with open(fname, 'rb') as fh: outstr += fh.read() elif fname.endswith('.h5'): @@ -56,6 +56,6 @@ class PlotTestHarness(TestHarness): def test_plot(): - harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', + harness = PlotTestHarness(('plot_1.png', 'plot_2.png', 'plot_3.png', 'plot_4.h5')) harness.main() diff --git a/tests/regression_tests/plot_overlaps/results_true.dat b/tests/regression_tests/plot_overlaps/results_true.dat index 90ec8924a..cb04daaa4 100644 --- a/tests/regression_tests/plot_overlaps/results_true.dat +++ b/tests/regression_tests/plot_overlaps/results_true.dat @@ -1 +1 @@ -566103831cb8273b0578565c39d30e479664e2b1783d877b45b42cf4f3af0b01671b6db423114b09a74bbe1ddf51a7db565ff2118d6d1ee987b52318773b719a \ No newline at end of file +926065ceb2a9b8292fe6270317c38c4373473cfea19d2a8392a32e5ece8e314c04b9f032921d987bd195ae4b6f674d359b0e38302e6ae4c93b4ac9573a384ac6 \ No newline at end of file diff --git a/tests/regression_tests/plot_overlaps/test.py b/tests/regression_tests/plot_overlaps/test.py index 62236081b..cf23abc27 100644 --- a/tests/regression_tests/plot_overlaps/test.py +++ b/tests/regression_tests/plot_overlaps/test.py @@ -19,7 +19,7 @@ class PlotTestHarness(TestHarness): openmc.plot_geometry(openmc_exec=config['exe']) def _test_output_created(self): - """Make sure *.ppm has been created.""" + """Make sure *.png has been created.""" for fname in self._plot_names: assert os.path.exists(fname), 'Plot output file does not exist.' @@ -34,8 +34,8 @@ class PlotTestHarness(TestHarness): outstr = bytes() for fname in self._plot_names: - if fname.endswith('.ppm'): - # Add PPM output to results + if fname.endswith('.png'): + # Add PNG output to results with open(fname, 'rb') as fh: outstr += fh.read() elif fname.endswith('.h5'): @@ -56,6 +56,6 @@ class PlotTestHarness(TestHarness): def test_plot_overlap(): - harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', + harness = PlotTestHarness(('plot_1.png', 'plot_2.png', 'plot_3.png', 'plot_4.h5')) harness.main() diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py index 3b7a6bef2..c4fa20c80 100644 --- a/tests/regression_tests/plot_voxel/test.py +++ b/tests/regression_tests/plot_voxel/test.py @@ -25,7 +25,7 @@ class PlotVoxelTestHarness(TestHarness): glob.glob('plot_4.h5')) def _test_output_created(self): - """Make sure *.ppm has been created.""" + """Make sure plots have been created.""" for fname in self._plot_names: assert os.path.exists(fname), 'Plot output file does not exist.' diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index de81a4179..eb60821c6 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -76,4 +76,27 @@ 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + + + -2.0 0.0 2.0 0.2 0.3 0.2 + + + + + + + + + + + + + 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 66b610b3a..1815324f6 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.033600E-01 1.676108E-03 +2.980096E-01 9.632798E-04 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 91316acc6..006358ca7 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -55,18 +55,20 @@ class SourceTestHarness(PyAPITestHarness): energy1 = openmc.stats.Maxwell(1.2895e6) energy2 = openmc.stats.Watt(0.988e6, 2.249e-6) energy3 = openmc.stats.Tabular(E, p, interpolation='histogram') + energy4 = openmc.stats.Mixture([1, 2, 3], [energy1, energy2, energy3]) source1 = openmc.Source(spatial1, angle1, energy1, strength=0.5) source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3) source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1) source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1) source5 = openmc.Source(spatial5, angle3, energy3, strength=0.1) + source6 = openmc.Source(spatial5, angle3, energy4, strength=0.1) settings = openmc.Settings() settings.batches = 10 settings.inactive = 5 settings.particles = 1000 - settings.source = [source1, source2, source3, source4, source5] + settings.source = [source1, source2, source3, source4, source5, source6] settings.export_to_xml() diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index cb4560d51..06b8e6bed 100644 --- a/tests/unit_tests/test_data_decay.py +++ b/tests/unit_tests/test_data_decay.py @@ -31,6 +31,18 @@ def u235_yields(): return openmc.data.FissionProductYields.from_endf(filename) +def test_get_decay_modes(): + assert openmc.data.get_decay_modes(1.0) == ['beta-'] + assert openmc.data.get_decay_modes(6.0) == ['sf'] + assert openmc.data.get_decay_modes(10.0) == ['unknown'] + + assert openmc.data.get_decay_modes(1.5) == ['beta-', 'n'] + assert openmc.data.get_decay_modes(1.4) == ['beta-', 'alpha'] + assert openmc.data.get_decay_modes(1.55) == ['beta-', 'n', 'n'] + assert openmc.data.get_decay_modes(1.555) == ['beta-', 'n', 'n', 'n'] + assert openmc.data.get_decay_modes(2.4) == ['ec/beta+', 'alpha'] + + def test_nb90_halflife(nb90): ufloat_close(nb90.half_life, ufloat(52560.0, 180.0)) ufloat_close(nb90.decay_constant, log(2.)/nb90.half_life) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index fa90b35f2..9769a49e3 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -7,7 +7,8 @@ import os from pathlib import Path import numpy as np -from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram, pool +from openmc.mpi import comm +from openmc.deplete import Chain, reaction_rates, nuclide, cram, pool import pytest from tests import cdtemp diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 3195c5976..a08d8738c 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -14,11 +14,11 @@ import numpy as np from uncertainties import ufloat import pytest +from openmc.mpi import comm from openmc.deplete import ( - ReactionRates, Results, ResultsList, comm, OperatorResult, - PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator, - EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, - cram) + ReactionRates, Results, ResultsList, OperatorResult, PredictorIntegrator, + CECMIntegrator, CF4Integrator, CELIIntegrator, EPCRK4Integrator, + LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram) from tests import dummy_operator diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 41613832f..ccff33ac3 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -263,6 +263,7 @@ def test_settings(lib_init): assert settings.generations_per_batch == 1 assert settings.particles == 100 assert settings.seed == 1 + assert settings.event_based is False settings.seed = 11 @@ -713,7 +714,6 @@ def test_trigger_set_n_batches(uo2_trigger_model, mpi_intracomm): def test_cell_translation(pincell_model_w_univ, mpi_intracomm): openmc.lib.finalize() openmc.lib.init(intracomm=mpi_intracomm) - openmc.lib.simulation_init() # Cell 1 is filled with a material so it has a translation, but we can't # set it. cell = openmc.lib.cells[1] @@ -727,9 +727,12 @@ def test_cell_translation(pincell_model_w_univ, mpi_intracomm): # This time we *can* set it cell.translation = (1., 0., -1.) assert cell.translation == pytest.approx([1., 0., -1.]) + openmc.lib.finalize() -def test_cell_rotation(pincell_model_w_univ): +def test_cell_rotation(pincell_model_w_univ, mpi_intracomm): + openmc.lib.finalize() + openmc.lib.init(intracomm=mpi_intracomm) # Cell 1 is filled with a material so we cannot rotate it, but we can get # its rotation matrix (which will be the identity matrix) cell = openmc.lib.cells[1] @@ -742,3 +745,4 @@ def test_cell_rotation(pincell_model_w_univ): assert cell.rotation == pytest.approx([0., 0., 0.]) cell.rotation = (180., 0., 0.) assert cell.rotation == pytest.approx([180., 0., 0.]) + openmc.lib.finalize() diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 4b658d5ab..d22485b3c 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -1,29 +1,281 @@ +from math import pi +from pathlib import Path +from shutil import which + +import numpy as np import pytest + import openmc import openmc.lib +@pytest.fixture(scope='function') +def pin_model_attributes(): + uo2 = openmc.Material(material_id=1, name='UO2') + uo2.set_density('g/cm3', 10.29769) + uo2.add_element('U', 1., enrichment=2.4) + uo2.add_element('O', 2.) + uo2.depletable = True + + zirc = openmc.Material(material_id=2, name='Zirc') + zirc.set_density('g/cm3', 6.55) + zirc.add_element('Zr', 1.) + zirc.depletable = False + + borated_water = openmc.Material(material_id=3, name='Borated water') + borated_water.set_density('g/cm3', 0.740582) + borated_water.add_element('B', 4.0e-5) + borated_water.add_element('H', 5.0e-2) + borated_water.add_element('O', 2.4e-2) + borated_water.add_s_alpha_beta('c_H_in_H2O') + borated_water.depletable = False + + mats = openmc.Materials([uo2, zirc, borated_water]) + + pitch = 1.25984 + fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR') + clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR') + box = openmc.model.rectangular_prism(pitch, pitch, + boundary_type='reflective') + + # Define cells + fuel_inf_cell = openmc.Cell(cell_id=1, name='inf fuel', fill=uo2) + fuel_inf_univ = openmc.Universe(universe_id=1, cells=[fuel_inf_cell]) + fuel = openmc.Cell(cell_id=2, name='fuel', + fill=fuel_inf_univ, region=-fuel_or) + clad = openmc.Cell(cell_id=3, fill=zirc, region=+fuel_or & -clad_or) + water = openmc.Cell(cell_id=4, fill=borated_water, region=+clad_or & box) + + # Define overall geometry + geom = openmc.Geometry([fuel, clad, water]) + uo2.volume = pi * fuel_or.r**2 + + settings = openmc.Settings() + settings.batches = 100 + settings.inactive = 10 + settings.particles = 1000 + + # Create a uniform spatial source distribution over fissionable zones + bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] + uniform_dist = openmc.stats.Box( + bounds[:3], bounds[3:], only_fissionable=True) + settings.source = openmc.source.Source(space=uniform_dist) + + entropy_mesh = openmc.RegularMesh() + entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] + entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] + entropy_mesh.dimension = [10, 10, 1] + settings.entropy_mesh = entropy_mesh + + tals = openmc.Tallies() + tal = openmc.Tally(tally_id=1, name='test') + tal.filters = [openmc.MaterialFilter(bins=[uo2])] + tal.scores = ['flux', 'fission'] + tals.append(tal) + + plot1 = openmc.Plot(plot_id=1) + plot1.origin = (0., 0., 0.) + plot1.width = (pitch, pitch) + plot1.pixels = (300, 300) + plot1.color_by = 'material' + plot1.filename = 'test' + plot2 = openmc.Plot(plot_id=2) + plot2.origin = (0., 0., 0.) + plot2.width = (pitch, pitch) + plot2.pixels = (300, 300) + plot2.color_by = 'cell' + plots = openmc.Plots((plot1, plot2)) + + chain = './test_chain.xml' + + chain_file_xml = """ + + + + + + 2.53000e-02 + + Xe136 + 1.0 + + + + +""" + operator_kwargs = {'chain_file': chain} + + return (mats, geom, settings, tals, plots, operator_kwargs, chain_file_xml) + + +def test_init(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, _, _ = \ + pin_model_attributes + + openmc.reset_auto_ids() + # Check blank initialization of a model + test_model = openmc.Model() + assert test_model.geometry.root_universe is None + assert len(test_model.materials) == 0 + ref_settings = openmc.Settings() + assert sorted(test_model.settings.__dict__.keys()) == \ + sorted(ref_settings.__dict__.keys()) + for ref_k, ref_v in ref_settings.__dict__.items(): + assert test_model.settings.__dict__[ref_k] == ref_v + assert len(test_model.tallies) == 0 + assert len(test_model.plots) == 0 + assert test_model._materials_by_id == {} + assert test_model._materials_by_name == {} + assert test_model._cells_by_id == {} + assert test_model._cells_by_name == {} + assert test_model.is_initialized is False + + # Now check proper init of an actual model. Assume no interference between + # parameters and so we can apply them all at once instead of testing one + # parameter initialization at a time + test_model = openmc.Model(geom, mats, settings, tals, plots) + assert test_model.geometry is geom + assert test_model.materials is mats + assert test_model.settings is settings + assert test_model.tallies is tals + assert test_model.plots is plots + assert test_model._materials_by_id == {1: mats[0], 2: mats[1], 3: mats[2]} + assert test_model._materials_by_name == { + 'UO2': {mats[0]}, 'Zirc': {mats[1]}, 'Borated water': {mats[2]}} + # The last cell is the one that contains the infinite fuel + assert test_model._cells_by_id == \ + {2: geom.root_universe.cells[2], 3: geom.root_universe.cells[3], + 4: geom.root_universe.cells[4], + 1: geom.root_universe.cells[2].fill.cells[1]} + # No cell name for 2 and 3, so we expect a blank name to be assigned to + # cell 3 due to overwriting + assert test_model._cells_by_name == { + 'fuel': {geom.root_universe.cells[2]}, + '': {geom.root_universe.cells[3], geom.root_universe.cells[4]}, + 'inf fuel': {geom.root_universe.cells[2].fill.cells[1]}} + assert test_model.is_initialized is False + + # Finally test the parameter type checking by passing bad types and + # obtaining the right exception types + def_params = [geom, mats, settings, tals, plots] + for i in range(len(def_params)): + args = def_params.copy() + # Try an integer, as that is a bad type for all arguments + args[i] = i + with pytest.raises(TypeError): + test_model = openmc.Model(*args) + + +def test_from_xml(run_in_tmpdir, pin_model_attributes): + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + + # This test will write the individual files to xml and then init that way + # and run the same sort of test as in test_init + mats.export_to_xml() + geom.export_to_xml() + settings.export_to_xml() + tals.export_to_xml() + plots.export_to_xml() + + # This from_xml method cannot load chain and fission_q + test_model = openmc.Model.from_xml() + assert test_model.geometry.root_universe.cells.keys() == \ + geom.root_universe.cells.keys() + assert [c.fill.name for c in + test_model.geometry.root_universe.cells.values()] == \ + [c.fill.name for c in geom.root_universe.cells.values()] + assert [mat.name for mat in test_model.materials] == \ + [mat.name for mat in mats] + # We will assume the attributes of settings that are custom objects are + # OK if the others are so we dotn need to implement explicit comparisons + no_test = ['_source', '_entropy_mesh'] + assert sorted(k for k in test_model.settings.__dict__.keys() + if k not in no_test) == \ + sorted(k for k in settings.__dict__.keys() if k not in no_test) + keys = sorted(k for k in settings.__dict__.keys() if k not in no_test) + for ref_k in keys: + assert test_model.settings.__dict__[ref_k] == settings.__dict__[ref_k] + assert len(test_model.tallies) == 0 + assert len(test_model.plots) == 0 + assert test_model._materials_by_id == \ + {1: test_model.materials[0], 2: test_model.materials[1], + 3: test_model.materials[2]} + assert test_model._materials_by_name == { + 'UO2': {test_model.materials[0]}, 'Zirc': {test_model.materials[1]}, + 'Borated water': {test_model.materials[2]}} + assert test_model._cells_by_id == { + 2: test_model.geometry.root_universe.cells[2], + 3: test_model.geometry.root_universe.cells[3], + 4: test_model.geometry.root_universe.cells[4], + 1: test_model.geometry.root_universe.cells[2].fill.cells[1]} + # No cell name for 2 and 3, so we expect a blank name to be assigned to + # cell 3 due to overwriting + assert test_model._cells_by_name == { + 'fuel': {test_model.geometry.root_universe.cells[2]}, + '': {test_model.geometry.root_universe.cells[3], + test_model.geometry.root_universe.cells[4]}, + 'inf fuel': {test_model.geometry.root_universe.cells[2].fill.cells[1]}} + assert test_model.is_initialized is False + + +def test_init_finalize_lib(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + # We are going to init and then make sure data is loaded + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + test_model = openmc.Model(geom, mats, settings, tals, plots) + test_model.init_lib(output=False, intracomm=mpi_intracomm) + + # First check that the API is advertised as initialized + assert openmc.lib.is_initialized is True + assert test_model.is_initialized is True + # Now make sure it actually is initialized by making a call to the lib + c_mat = openmc.lib.find_material((0.6, 0., 0.)) + # This should be Borated water + assert c_mat.name == 'Borated water' + assert c_mat.id == 3 + + # Ok, now lets test that we can clear the data and check that it is cleared + test_model.finalize_lib() + + # First check that the API is advertised as initialized + assert openmc.lib.is_initialized is False + assert test_model.is_initialized is False + # Note we cant actually test that a sys call fails because we should get a + # seg fault + + def test_import_properties(run_in_tmpdir, mpi_intracomm): """Test importing properties on the Model class """ # Create PWR pin cell model and write XML files openmc.reset_auto_ids() model = openmc.examples.pwr_pin_cell() - model.export_to_xml() + model.init_lib(output=False, intracomm=mpi_intracomm) # Change fuel temperature and density and export properties - openmc.lib.init(intracomm=mpi_intracomm) cell = openmc.lib.cells[1] cell.set_temperature(600.0) cell.fill.set_density(5.0, 'g/cm3') - openmc.lib.export_properties() + openmc.lib.export_properties(output=False) + + # Import properties to existing model + model.import_properties("properties.h5") + + # Check to see that values are assigned to the C and python representations + # First python + cell = model.geometry.get_all_cells()[1] + assert cell.temperature == [600.0] + assert cell.fill.get_mass_density() == pytest.approx(5.0) + # Now C + assert openmc.lib.cells[1].get_temperature() == 600. + assert openmc.lib.materials[1].get_density('g/cm3') == pytest.approx(5.0) + + # Clear the C API openmc.lib.finalize() - # Import properties to existing model and re-export to new directory - model.import_properties("properties.h5") + # Verify the attributes survived by exporting to XML and re-creating model.export_to_xml("with_properties") - # Load model with properties and confirm temperature/density has been changed + # Load model with properties and confirm temperature/density changed model_with_properties = openmc.Model.from_xml( 'with_properties/geometry.xml', 'with_properties/materials.xml', @@ -32,3 +284,251 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm): cell = model_with_properties.geometry.get_all_cells()[1] assert cell.temperature == [600.0] assert cell.fill.get_mass_density() == pytest.approx(5.0) + + +def test_run(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + test_model = openmc.Model(geom, mats, settings, tals, plots) + + # This case will run by getting the k-eff and tallies for command-line and + # C API execution modes and ensuring they give the same result. + sp_path = test_model.run(output=False) + with openmc.StatePoint(sp_path) as sp: + cli_keff = sp.k_combined + cli_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0] + cli_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0] + + test_model.init_lib(output=False, intracomm=mpi_intracomm) + sp_path = test_model.run(output=False) + with openmc.StatePoint(sp_path) as sp: + lib_keff = sp.k_combined + lib_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0] + lib_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0] + + # and lets compare results + assert lib_keff.n == pytest.approx(cli_keff.n, abs=1e-13) + assert lib_flux == pytest.approx(cli_flux, abs=1e-13) + assert lib_fiss == pytest.approx(cli_fiss, abs=1e-13) + + # Now we should make sure that the flags for items which should be handled + # by init are properly set + with pytest.raises(ValueError): + test_model.run(threads=1) + with pytest.raises(ValueError): + test_model.run(geometry_debug=True) + with pytest.raises(ValueError): + test_model.run(restart_file='1.h5') + with pytest.raises(ValueError): + test_model.run(tracks=True) + + test_model.finalize_lib() + + +def test_plots(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + test_model = openmc.Model(geom, mats, settings, tals, plots) + + # This test cannot check the correctness of the plot, but it can + # check that a plot was made and that the expected png files are there + + # We will run the test twice, the first time without C API, the second with + for i in range(2): + if i == 1: + test_model.init_lib(output=False, intracomm=mpi_intracomm) + test_model.plot_geometry(output=False) + + # Now look for the files + for fname in ('test.png', 'plot_2.png'): + test_file = Path(fname) + assert test_file.exists() + test_file.unlink() + + test_model.finalize_lib() + + +def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + test_model = openmc.Model(geom, mats, settings, tals, plots) + + test_model.init_lib(output=False, intracomm=mpi_intracomm) + + # Now we can call rotate_cells, translate_cells, update_densities, + # and update_cell_temperatures and make sure the changes have taken hold. + # For each we will first try bad inputs to make sure we get the right + # errors and then we do a good one which calls the material by name and + # then id to make sure it worked + + # The rotate_cells and translate_cells will work on the cell named fill, as + # it is filled with a universe and thus the operation will be valid + + # First rotate_cells + with pytest.raises(TypeError): + # Make sure it tells us we have a bad names_or_ids type + test_model.rotate_cells(None, (0, 0, 90)) + with pytest.raises(TypeError): + test_model.rotate_cells([None], (0, 0, 90)) + with pytest.raises(openmc.exceptions.InvalidIDError): + # Make sure it tells us we had a bad id + test_model.rotate_cells([7200], (0, 0, 90)) + with pytest.raises(openmc.exceptions.InvalidIDError): + # Make sure it tells us we had a bad id + test_model.rotate_cells(['bad_name'], (0, 0, 90)) + # Now a good one + assert np.all(openmc.lib.cells[2].rotation == (0., 0., 0.)) + test_model.rotate_cells([2], (0, 0, 90)) + assert np.all(openmc.lib.cells[2].rotation == (0., 0., 90.)) + + # And same thing by name + test_model.rotate_cells(['fuel'], (0, 0, 180)) + + # Now translate_cells. We dont need to re-check the TypeErrors/bad ids, + # because the other functions use the same hidden method as rotate_cells + assert np.all(openmc.lib.cells[2].translation == (0., 0., 0.)) + test_model.translate_cells([2], (0, 0, 10)) + assert np.all(openmc.lib.cells[2].translation == (0., 0., 10.)) + + # Now lets do the density updates. + # Check initial conditions + assert openmc.lib.materials[1].get_density('atom/b-cm') == \ + pytest.approx(0.06891296988603757, abs=1e-13) + mat_a_dens = np.sum( + [v[1] for v in test_model.materials[0]. + get_nuclide_atom_densities().values()]) + assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8) + # Change the density + test_model.update_densities(['UO2'], 2.) + assert openmc.lib.materials[1].get_density('atom/b-cm') == \ + pytest.approx(2., abs=1e-13) + mat_a_dens = np.sum( + [v[1] for v in test_model.materials[0]. + get_nuclide_atom_densities().values()]) + assert mat_a_dens == pytest.approx(2., abs=1e-8) + + # Now lets do the cell temperature updates. + # Check initial conditions + assert test_model._cells_by_id == \ + {2: geom.root_universe.cells[2], 3: geom.root_universe.cells[3], + 4: geom.root_universe.cells[4], + 1: geom.root_universe.cells[2].fill.cells[1]} + assert openmc.lib.cells[3].get_temperature() == \ + pytest.approx(293.6, abs=1e-13) + assert test_model.geometry.root_universe.cells[3].temperature is None + # Change the temperature + test_model.update_cell_temperatures([3], 600.) + assert openmc.lib.cells[3].get_temperature() == \ + pytest.approx(600., abs=1e-13) + assert test_model.geometry.root_universe.cells[3].temperature == \ + pytest.approx(600., abs=1e-13) + + # And finally material volume + assert openmc.lib.materials[1].volume == \ + pytest.approx(0.4831931368640985, abs=1e-13) + # The temperature on the material will be None because its just the default + assert test_model.materials[0].volume == \ + pytest.approx(0.4831931368640985, abs=1e-13) + # Change the temperature + test_model.update_material_volumes(['UO2'], 2.) + assert openmc.lib.materials[1].volume == pytest.approx(2., abs=1e-13) + assert test_model.materials[0].volume == pytest.approx(2., abs=1e-13) + + test_model.finalize_lib() + + +def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, op_kwargs, chain_file_xml = \ + pin_model_attributes + with open('test_chain.xml', 'w') as f: + f.write(chain_file_xml) + test_model = openmc.Model(geom, mats, settings, tals, plots) + + initial_mat = mats[0].clone() + initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1] + + # Note that the chain file includes only U-235 fission to a stable Xe136 w/ + # a yield of 100%. Thus all the U235 we lose becomes Xe136 + + # In this test we first run without pre-initializing the shared library + # data and then compare. Then we repeat with the C API already initialized + # and make sure we get the same answer + test_model.deplete([1e6], 'predictor', final_step=False, + operator_kwargs=op_kwargs, + power=1., output=False) + # Get the new Xe136 and U235 atom densities + after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] + after_u = mats[0].get_nuclide_atom_densities()['U235'][1] + assert after_xe + after_u == pytest.approx(initial_u, abs=1e-15) + assert test_model.is_initialized is False + + # Reset the initial material densities + mats[0].nuclides.clear() + densities = initial_mat.get_nuclide_atom_densities() + tot_density = 0. + for nuc, density in densities.values(): + mats[0].add_nuclide(nuc, density) + tot_density += density + mats[0].set_density('atom/b-cm', tot_density) + + # Now we can re-run with the pre-initialized API + test_model.init_lib(output=False, intracomm=mpi_intracomm) + test_model.deplete([1e6], 'predictor', final_step=False, + operator_kwargs=op_kwargs, + power=1., output=False) + # Get the new Xe136 and U235 atom densities + after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] + after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1] + assert after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15) + assert test_model.is_initialized is True + + # And end by comparing to the previous case + assert after_xe == pytest.approx(after_lib_xe, abs=1e-15) + assert after_u == pytest.approx(after_lib_u, abs=1e-15) + + test_model.finalize_lib() + + +def test_calc_volumes(run_in_tmpdir, pin_model_attributes, mpi_intracomm): + mats, geom, settings, tals, plots, _, _ = pin_model_attributes + + test_model = openmc.Model(geom, mats, settings, tals, plots) + + # With no vol calcs, it should fail + with pytest.raises(ValueError): + test_model.calculate_volumes(output=False) + + # Add a cell and mat volume calc + material_vol_calc = openmc.VolumeCalculation( + [mats[2]], samples=1000, lower_left=(-.63, -.63, -100.), + upper_right=(.63, .63, 100.)) + cell_vol_calc = openmc.VolumeCalculation( + [geom.root_universe.cells[3]], samples=1000, + lower_left=(-.63, -.63, -100.), upper_right=(.63, .63, 100.)) + test_model.settings.volume_calculations = \ + [material_vol_calc, cell_vol_calc] + + # Now lets compute the volumes and check to see if it was applied + # First lets do without using the C-API + # Make sure the volumes are unassigned first + assert mats[2].volume is None + assert geom.root_universe.cells[3].volume is None + test_model.calculate_volumes(output=False, apply_volumes=True) + + # Now let's test that we have volumes assigned; we arent checking the + # value, just that the value was changed + assert mats[2].volume > 0. + assert geom.root_universe.cells[3].volume > 0. + + # Now reset the values + mats[2].volume = None + geom.root_universe.cells[3].volume = None + + # And do again with an initialized library + for file in ['volume_1.h5', 'volume_2.h5']: + file = Path(file) + file.unlink() + test_model.init_lib(output=False, intracomm=mpi_intracomm) + test_model.calculate_volumes(output=False, apply_volumes=True) + assert mats[2].volume > 0. + assert geom.root_universe.cells[3].volume > 0. + assert openmc.lib.materials[3].volume == mats[2].volume + + test_model.finalize_lib() diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index e15958478..9e6237ba1 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -101,8 +101,12 @@ def test_mixture(): assert mix.distribution == [d1, d2] assert len(mix) == 4 - with pytest.raises(NotImplementedError): - mix.to_xml_element('distribution') + elem = mix.to_xml_element('distribution') + + d = openmc.stats.Mixture.from_xml_element(elem) + assert d.probability == p + assert d.distribution == [d1, d2] + assert len(d) == 4 def test_polar_azimuthal(): diff --git a/vendor/fmt b/vendor/fmt index 65ac626c5..d141cdbeb 160000 --- a/vendor/fmt +++ b/vendor/fmt @@ -1 +1 @@ -Subproject commit 65ac626c5856f5aad1f1542e79407a6714357043 +Subproject commit d141cdbeb0fb422a3fb7173b285fd38e0d1772dc