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

Merge remote-tracking branch 'upstream/develop' into cmfd-capi

Browse source
This commit is contained in:
Shikhar Kumar 2018-11-23 19:00:50 -05:00
commit 0dfaaea540
142 changed files with 3042 additions and 2519 deletions
Download patch file Download diff file Expand all files Collapse all files

6
.travis.yml
View file

@ -15,11 +15,12 @@ addons:
- libhdf5-mpich-dev
- libblas-dev
- liblapack-dev
config:
retries: true
cache:
directories:
- $HOME/nndc_hdf5
- $HOME/endf-b-vii.1
- $HOME/WMP_Library
env:
global:
- FC=gfortran
@ -28,7 +29,6 @@ env:
- OMP_NUM_THREADS=2
- OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
- OPENMC_MULTIPOLE_LIBRARY=$HOME/WMP_Library
- LD_LIBRARY_PATH=$HOME/MOAB/lib:$HOME/DAGMC/lib
- PATH=$PATH:$HOME/NJOY2016/build
- DISPLAY=:99.0
@ -39,7 +39,7 @@ env:
- OMP=n MPI=y PHDF5=n
- OMP=n MPI=y PHDF5=y
- OMP=n MPI=y PHDF5=y DAGMC=y
- OMP=y MPI=y PHDF5=y DAGMC=y
- OMP=y MPI=y PHDF5=y DAGMC=y
notifications:
webhooks: https://coveralls.io/webhook?repo_token=$COVERALLS_REPO_TOKEN
install:

3
CMakeLists.txt
View file

@ -385,6 +385,7 @@ add_library(libopenmc SHARED
src/cell.cpp
src/cmfd_execute.cpp
src/cmfd_solver.cpp
src/cross_sections.cpp
src/distribution.cpp
src/distribution_angle.cpp
src/distribution_energy.cpp
@ -426,7 +427,7 @@ add_library(libopenmc SHARED
src/simulation.cpp
src/source.cpp
src/state_point.cpp
src/string_functions.cpp
src/string_utils.cpp
src/summary.cpp
src/surface.cpp
src/tallies/filter.cpp

3
Dockerfile
View file

@ -5,7 +5,6 @@ ENV FC=/usr/bin/mpif90 CC=/usr/bin/mpicc CXX=/usr/bin/mpicxx \
PATH=/opt/openmc/bin:/opt/NJOY2016/build:$PATH \
LD_LIBRARY_PATH=/opt/openmc/lib:$LD_LIBRARY_PATH \
OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \
OPENMC_MULTIPOLE_LIBRARY=/root/WMP_Library \
OPENMC_ENDF_DATA=/root/endf-b-vii.1
# Install dependencies from Debian package manager
@ -36,4 +35,4 @@ RUN git clone https://github.com/openmc-dev/openmc.git /opt/openmc && \
cd .. && pip install -e .[test]
# Download cross sections (NNDC and WMP) and ENDF data needed by test suite
RUN ./opt/openmc/tools/ci/download-xs.sh
RUN ./opt/openmc/tools/ci/download-xs.sh

10
docs/source/capi/index.rst
View file

@ -333,12 +333,14 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_material_set_density(int32_t index, double density)
.. c:function:: int openmc_material_set_density(int32_t index, double density, const char* units)
Set the density of a material.
:param int32_t index: Index in the materials array
:param double density: Density of the material in atom/b-cm
:param double density: Density of the material
:param units: Units for density
:type units: const char*
:return: Return status (negative if an error occurs)
:rtype: int
@ -463,13 +465,15 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_statepoint_write(const char filename[])
.. c:function:: int openmc_statepoint_write(const char filename[], const bool* write_source)
Write a statepoint file
:param filename: Name of file to create. If a null pointer is passed, a
filename is assigned automatically.
:type filename: const char[]
:param write_source: Whether to include the source bank
:type write_source: const bool*
:return: Return status (negative if an error occurs)
:rtype: int

20
docs/source/devguide/styleguide.rst
View file

@ -237,6 +237,26 @@ are written differently by convention (e.g., ``E`` for energy). Data members of
classes (but not structs) additionally have trailing underscores (e.g.,
``a_class_member_``).
The following conventions are used for variables with short names:
- ``d`` stands for "distance"
- ``E`` stands for "energy"
- ``p`` stands for "particle"
- ``r`` stands for "position"
- ``rx`` stands for "reaction"
- ``u`` stands for "direction"
- ``xs`` stands for "cross section"
All classes and non-member functions should be declared within the ``openmc``
namespace. Global variables must be declared in a namespace nested within the
``openmc`` namespace. The following sub-namespaces are in use:
- ``openmc::data``: Fundamental nuclear data (cross sections, multigroup data,
decay constants, etc.)
- ``openmc::model``: Variables related to geometry, materials, and tallies
- ``openmc::settings``: Global settings / options
- ``openmc::simulation``: Variables used only during a simulation
Accessors and mutators (get and set functions) may be named like
variables. These often correspond to actual member variables, but this is not
required. For example, ``int count()`` and ``void set_count(int count)``.

10
docs/source/devguide/tests.rst
View file

@ -30,11 +30,11 @@ or using pip (recommended)::
pip install -e .[test]
It is also assumed that you have cross section data available that is pointed to
by the :envvar:`OPENMC_CROSS_SECTIONS` and :envvar:`OPENMC_MULTIPOLE_LIBRARY`
environment variables. Furthermore, to run unit tests for the :mod:`openmc.data`
module, it is necessary to have ENDF/B-VII.1 data available and pointed to by
the :envvar:`OPENMC_ENDF_DATA` environment variable. All data sources can be
obtained using the ``tools/ci/travis-before-script.sh`` script.
by the :envvar:`OPENMC_CROSS_SECTIONS` environment variables. Furthermore, to
run unit tests for the :mod:`openmc.data` module, it is necessary to have
ENDF/B-VII.1 data available and pointed to by the :envvar:`OPENMC_ENDF_DATA`
environment variable. All data sources can be obtained using the
``tools/ci/travis-before-script.sh`` script.
To execute the test suite, go to the ``tests/`` directory and run::

10
docs/source/io_formats/cross_sections.rst
View file

@ -25,9 +25,10 @@ node. For example,
``<library>`` Element
---------------------
The ``<library>`` element indicates where an HDF5 cross section file is located,
whether it contains incident neutron or thermal scattering data, and what
materials are listed within. It has the following attributes:
The ``<library>`` element indicates where an HDF5 data file is located, whether
it contains incident neutron, incident photon, thermal scattering, or windowed
multipole data, and what materials are listed within. It has the following
attributes:
:materials:
@ -48,4 +49,5 @@ materials are listed within. It has the following attributes:
directory containing the ``cross_sections.xml`` file.
:type:
The type of data contained in the file, either 'neutron' or 'thermal'.
The type of data contained in the file. Accepted values are 'neutron',
'thermal', 'photon', and 'wmp'.

15
docs/source/io_formats/materials.rst
View file

@ -18,21 +18,6 @@ path to the XML cross section listing when in continuous-energy mode, and the
:envvar:`OPENMC_MG_CROSS_SECTIONS` environment variable will be used in
multi-group mode.
.. _multipole_library:
-------------------------------
``<multipole_library>`` Element
-------------------------------
The ``<multipole_library>`` element indicates the directory containing a
windowed multipole library. If a windowed multipole library is available,
OpenMC can use it for on-the-fly Doppler-broadening of resolved resonance range
cross sections. If this element is absent from the settings.xml file, the
:envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used.
.. note:: The <temperature_multipole> element must also be set to "true" for
windowed multipole functionality.
.. _material:
----------------------

2
docs/source/pythonapi/index.rst
View file

@ -22,7 +22,7 @@ there are many substantial benefits to using the Python API, including:
- Ability to plot individual universes as geometry is being created
- A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`)
- Random sphere packing for generating TRISO particle locations
(:func:`openmc.model.pack_trisos`)
(:func:`openmc.model.pack_spheres`)
- Ability to create materials based on natural elements or uranium enrichment
For those new to Python, there are many good tutorials available online. We

2
docs/source/pythonapi/model.rst
View file

@ -37,7 +37,7 @@ Functions
:template: myfunction.rst
openmc.model.create_triso_lattice
openmc.model.pack_trisos
openmc.model.pack_spheres
Model Container
---------------

17
docs/source/usersguide/cross_sections.rst
View file

@ -41,11 +41,6 @@ following environment variables are used:
user has not specified :attr:`Materials.cross_sections` (equivalently, the
:ref:`cross_sections` in :ref:`materials.xml <io_materials>`).
:envvar:`OPENMC_MULTIPOLE_LIBRARY`
Indicates the path to a directory containing windowed multipole data if the
user has not specified :attr:`Materials.multipole_library` (equivalently, the
:ref:`multipole_library` in :ref:`materials.xml <io_materials>`)
:envvar:`OPENMC_MG_CROSS_SECTIONS`
Indicates the path to the an :ref:`HDF5 file <io_mgxs_library>` that contains
multi-group cross sections if the user has not specified
@ -306,12 +301,12 @@ Windowed Multipole Data
-----------------------
OpenMC is capable of using windowed multipole data for on-the-fly Doppler
broadening. While such data is not yet available for all nuclides, an
experimental multipole library is available that contains data for 70
nuclides. To obtain this library, you can run :ref:`scripts_multipole` which
will download and extract it into a ``wmp`` directory. Once the library has been
downloaded, set the :envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable (or
the :attr:`Materials.multipole_library` attribute) to the ``wmp`` directory.
broadening. A comprehensive multipole data library containing all nuclides in
ENDF/B-VII.1 is available on `GitHub
<https://github.com/mit-crpg/WMP_Library>`_. To obtain this library, download
and unpack an archive (.zip or .tag.gz) from GitHub. Once unpacked, you can use
the :class:`openmc.data.DataLibrary` class to register the .h5 files as
described in :ref:`create_xs_library`.
--------------------------
Multi-Group Cross Sections

11
docs/source/usersguide/scripts.rst
View file

@ -149,17 +149,6 @@ the following optional arguments:
and processing the data may require as much as 40 GB of additional
free disk space.
.. _scripts_multipole:
-----------------------------
``openmc-get-multipole-data``
-----------------------------
This script downloads and extracts windowed multipole data based on
ENDF/B-VII.1. It has the following optional arguments:
-b, --batch Suppress standard in
.. _scripts_nndc:
------------------------

524
examples/jupyter/nuclear-data.ipynb
View file

File diff suppressed because one or more lines are too long

115
examples/jupyter/triso.ipynb
View file

File diff suppressed because one or more lines are too long

9
include/openmc/capi.h
View file

@ -60,7 +60,7 @@ extern "C" {
int openmc_material_get_id(int32_t index, int32_t* id);
int openmc_material_get_fissionable(int32_t index, bool* fissionable);
int openmc_material_get_volume(int32_t index, double* volume);
int openmc_material_set_density(int32_t index, double density);
int openmc_material_set_density(int32_t index, double density, const char* units);
int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density);
int openmc_material_set_id(int32_t index, int32_t id);
int openmc_material_set_volume(int32_t index, double volume);
@ -135,17 +135,10 @@ extern "C" {
// Global variables
extern char openmc_err_msg[256];
extern int32_t n_cells;
extern int32_t n_lattices;
extern int32_t n_materials;
extern int n_nuclides;
extern int32_t n_plots;
extern int32_t n_realizations;
extern int32_t n_sab_tables;
extern int32_t n_sources;
extern int32_t n_surfaces;
extern int32_t n_tallies;
extern int32_t n_universes;
// Variables that are shared by necessity (can be removed from public header
// later)

17
include/openmc/cell.h
View file

@ -39,16 +39,21 @@ constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
// Global variables
//==============================================================================
class Cell;
class Universe;
namespace model {
extern "C" int32_t n_cells;
class Cell;
extern std::vector<Cell*> cells;
extern std::unordered_map<int32_t, int32_t> cell_map;
class Universe;
extern std::vector<Universe*> universes;
extern std::unordered_map<int32_t, int32_t> universe_map;
} // namespace model
//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================
@ -145,13 +150,13 @@ public:
virtual ~Cell() {}
};
class CSGCell : public Cell
{
public:
CSGCell();
explicit CSGCell(pugi::xml_node cell_node);
bool
@ -163,7 +168,7 @@ public:
void to_hdf5(hid_t group_id) const;
protected:
bool contains_simple(Position r, Direction u, int32_t on_surface) const;
bool contains_complex(Position r, Direction u, int32_t on_surface) const;
@ -183,6 +188,6 @@ public:
};
#endif
} // namespace openmc
#endif // OPENMC_CELL_H

2
include/openmc/constants.h
View file

@ -25,7 +25,7 @@ constexpr int VERSION_RELEASE {0};
constexpr std::array<int, 3> VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE};
// HDF5 data format
constexpr int HDF5_VERSION[] {1, 0};
constexpr int HDF5_VERSION[] {2, 0};
// Version numbers for binary files
constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};

67
include/openmc/cross_sections.h Normal file
View file

@ -0,0 +1,67 @@
#ifndef OPENMC_CROSS_SECTIONS_H
#define OPENMC_CROSS_SECTIONS_H
#include "pugixml.hpp"
#include <string>
#include <map>
#include <vector>
namespace openmc {
//==============================================================================
// Library class
//==============================================================================
class Library {
public:
// Types, enums
enum class Type {
neutron = 1, photon = 3, thermal = 2, multigroup = 4, wmp = 5
};
// Constructors
Library() { };
Library(pugi::xml_node node, const std::string& directory);
// Comparison operator (for using in map)
bool operator<(const Library& other) {
return path_ < other.path_;
}
// Data members
Type type_; //!< Type of data library
std::vector<std::string> materials_; //!< Materials contained in library
std::string path_; //!< File path to library
};
using LibraryKey = std::pair<Library::Type, std::string>;
//==============================================================================
// Global variable declarations
//==============================================================================
namespace data {
//!< Data libraries
extern std::vector<Library> libraries;
//! Maps (type, name) to index in libraries
extern std::map<LibraryKey, std::size_t> library_map;
} // namespace data
//==============================================================================
// Non-member functions
//==============================================================================
//! Read cross sections file (either XML or multigroup H5) and populate data
//! libraries
extern "C" void read_cross_sections_xml();
//! Read cross_sections.xml and populate data libraries
void read_ce_cross_sections_xml();
} // namespace openmc
#endif // OPENMC_CROSS_SECTIONS_H

14
include/openmc/dagmc.h
View file

@ -10,12 +10,24 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace model {
extern moab::DagMC* DAG;
} // namespace model
//==============================================================================
// Non-member functions
//==============================================================================
extern "C" void load_dagmc_geometry();
extern "C" void free_memory_dagmc();
}
} // namespace openmc
#endif // DAGMC

4
include/openmc/eigenvalue.h
View file

@ -18,6 +18,8 @@ namespace openmc {
// Global variables
//==============================================================================
namespace simulation {
extern double keff_generation; //!< Single-generation k on each processor
extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
extern std::vector<double> entropy; //!< Shannon entropy at each generation
@ -26,6 +28,8 @@ extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
extern "C" int64_t n_bank;
#pragma omp threadprivate(n_bank)
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================

6
include/openmc/geometry.h
View file

@ -13,10 +13,14 @@ namespace openmc {
// Global variables
//==============================================================================
extern "C" int openmc_root_universe;
namespace model {
extern "C" int root_universe;
extern std::vector<int64_t> overlap_check_count;
} // namespace model
//==============================================================================
//! Check for overlapping cells at a particle's position.
//==============================================================================

6
include/openmc/lattice.h
View file

@ -31,10 +31,14 @@ enum class LatticeType {
//==============================================================================
class Lattice;
extern std::vector<Lattice*> lattices;
namespace model {
extern std::vector<Lattice*> lattices;
extern std::unordered_map<int32_t, int32_t> lattice_map;
} // namespace model
//==============================================================================
//! \class Lattice
//! \brief Abstract type for ordered array of universes.

5
include/openmc/material.h
View file

@ -14,9 +14,14 @@ namespace openmc {
//==============================================================================
class Material;
namespace model {
extern std::vector<Material*> materials;
extern std::unordered_map<int32_t, int32_t> material_map;
} // namespace model
//==============================================================================
//! A substance with constituent nuclides and thermal scattering data
//==============================================================================

22
include/openmc/mesh.h
View file

@ -17,6 +17,19 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
class RegularMesh;
namespace model {
extern std::vector<std::unique_ptr<RegularMesh>> meshes;
extern std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
//==============================================================================
//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
//==============================================================================
@ -117,15 +130,6 @@ extern "C" void read_meshes(pugi::xml_node* root);
//! \param[in] group HDF5 group
extern "C" void meshes_to_hdf5(hid_t group);
//==============================================================================
// Global variables
//==============================================================================
extern std::vector<std::unique_ptr<RegularMesh>> meshes;
extern std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace openmc
#endif // OPENMC_MESH_H

4
include/openmc/mgxs_interface.h
View file

@ -15,6 +15,8 @@ namespace openmc {
// Global variables
//==============================================================================
namespace data {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
@ -22,6 +24,8 @@ extern std::vector<double> energy_bins;
extern std::vector<double> energy_bin_avg;
extern std::vector<double> rev_energy_bins;
} // namespace data
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================

4
include/openmc/nuclide.h
View file

@ -14,11 +14,15 @@ namespace openmc {
// Global variables
//==============================================================================
namespace data {
// Minimum/maximum transport energy for each particle type. Order corresponds to
// that of the ParticleType enum
extern std::array<double, 2> energy_min;
extern std::array<double, 2> energy_max;
} // namespace data
//===============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy

13
include/openmc/physics_mg.h
View file

@ -10,34 +10,29 @@
namespace openmc {
//TODO: Remove energy_bin_avg and material_xs parameters when they reside on
//TODO: Remove material_xs parameters when they reside on
// the C-side this should happen after materials, physics, input, and tallies
// are brought over
//! \brief samples particle behavior after a collision event.
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
//! \param material_xs The cross section cache for the current material
extern "C" void
collision_mg(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs);
collision_mg(Particle* p, const MaterialMacroXS* material_xs);
//! \brief samples a reaction type.
//!
//! Note that there is special logic when suvival biasing is turned on since
//! fission and disappearance are treated implicitly.
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
//! \param material_xs The cross section cache for the current material
void
sample_reaction(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs);
sample_reaction(Particle* p, const MaterialMacroXS* material_xs);
//! \brief Samples the scattering event
//! \param p Particle to operate on
//! \param energy_bin_avg Average energy within each energy bin
void
scatter(Particle* p, const double* energy_bin_avg);
scatter(Particle* p);
//! \brief Determines the average total, prompt and delayed neutrons produced
//! from fission and creates the appropriate bank sites.

18
include/openmc/plot.h
View file

@ -18,14 +18,14 @@ namespace openmc {
// Global variables
//===============================================================================
extern int PLOT_LEVEL_LOWEST; //!< lower bound on plot universe level
class Plot;
namespace model {
extern std::vector<Plot> plots; //!< Plot instance container
extern std::unordered_map<int, int> plot_map; //!< map of plot ids to index
extern "C" int32_t n_plots; //!< number of plots in openmc run
class Plot;
extern std::vector<Plot> plots; //!< Plot instance container
} // namespace model
//===============================================================================
// RGBColor holds color information for plotted objects
@ -124,8 +124,7 @@ void draw_mesh_lines(Plot pl, ImageData& data);
//! Write a ppm image to file using a plot object's image data
//! \param[in] plot object
//! \param[out] image data associated with the plot object
void output_ppm(Plot pl,
const ImageData& data);
void output_ppm(Plot pl, const ImageData& data);
//! Get the rgb color for a given particle position in a plot
//! \param[in] particle with position for current pixel
@ -141,7 +140,7 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id);
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to memory space of voxel data
void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
hid_t* dset, hid_t* memspace);
hid_t* dset, hid_t* memspace);
//! Write a section of the voxel data to hdf5
//! \param[in] voxel slice
@ -149,7 +148,8 @@ void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
//! \param[out] dataset pointer to voxesl data
//! \param[out] pointer to data to write
void voxel_write_slice(int x, hid_t dspace, hid_t dset,
hid_t memspace, void* buf);
hid_t memspace, void* buf);
//! Close voxel file entities
//! \param[in] data space to close
//! \param[in] dataset to close

4
include/openmc/settings.h
View file

@ -52,7 +52,6 @@ extern "C" bool dagmc; //!< indicator of DAGMC geometry
// Paths to various files
extern std::string path_cross_sections; //!< path to cross_sections.xml
extern std::string path_input; //!< directory where main .xml files resides
extern std::string path_multipole; //!< directory containing multipole files
extern std::string path_output; //!< directory where output files are written
extern std::string path_particle_restart; //!< path to a particle restart file
extern std::string path_source;
@ -62,7 +61,7 @@ extern std::string path_statepoint; //!< path to a statepoint file
extern "C" int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array
extern "C" int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array
extern "C" int32_t index_cmfd_mesh; //!< Index of CMFD mesh in global mesh array
extern "C" int32_t n_batches; //!< number of (inactive+active) batches
extern "C" int32_t n_inactive; //!< number of inactive batches
extern "C" int32_t gen_per_batch; //!< number of generations per batch
@ -77,6 +76,7 @@ extern "C" int n_max_batches; //!< Maximum number of batches
extern "C" int res_scat_method; //!< resonance upscattering method
extern "C" double res_scat_energy_min; //!< Min energy in [eV] for res. upscattering
extern "C" double res_scat_energy_max; //!< Max energy in [eV] for res. upscattering
extern std::vector<std::string> res_scat_nuclides; //!< Nuclides using res. upscattering treatment
extern "C" int run_mode; //!< Run mode (eigenvalue, fixed src, etc.)
extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written

18
include/openmc/source.h
View file

@ -16,6 +16,18 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
class SourceDistribution;
namespace model {
extern std::vector<SourceDistribution> external_sources;
} // namespace model
//==============================================================================
//! External source distribution
//==============================================================================
@ -40,12 +52,6 @@ private:
UPtrDist energy_; //!< Energy distribution
};
//==============================================================================
// Global variables
//==============================================================================
extern std::vector<SourceDistribution> external_sources;
//==============================================================================
// Functions
//==============================================================================

20
include/openmc/string_functions.h
View file

@ -1,20 +0,0 @@
//! \file string_functions.h
//! A collection of helper routines for C-strings and STL strings
#ifndef OPENMC_STRING_FUNCTIONS_H
#define OPENMC_STRING_FUNCTIONS_H
#include <string>
namespace openmc {
std::string& strtrim(std::string& s);
char* strtrim(char* c_str);
void to_lower(std::string& str);
int word_count(std::string const& str);
} // namespace openmc
#endif // STRING_FUNCTIONS_H

34
include/openmc/string_utils.h
View file

@ -1,42 +1,24 @@
#ifndef OPENMC_STRING_UTILS_H
#define OPENMC_STRING_UTILS_H
#include <algorithm>
#include <string>
#include <vector>
namespace openmc {
inline std::vector<std::string>
split(const std::string& in)
{
std::vector<std::string> out;
std::string& strtrim(std::string& s);
for (int i = 0; i < in.size(); ) {
// Increment i until we find a non-whitespace character.
if (std::isspace(in[i])) {
i++;
char* strtrim(char* c_str);
} else {
// Find the next whitespace character at j.
int j = i + 1;
while (j < in.size() && std::isspace(in[j]) == 0) {j++;}
void to_lower(std::string& str);
// Push-back everything between i and j.
out.push_back(in.substr(i, j-i));
i = j + 1; // j is whitespace so leapfrog to j+1
}
}
int word_count(std::string const& str);
return out;
}
std::vector<std::string> split(const std::string& in);
inline bool
ends_with(const std::string& value, const std::string& ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
bool ends_with(const std::string& value, const std::string& ending);
bool starts_with(const std::string& value, const std::string& beginning);
} // namespace openmc
#endif // OPENMC_STRING_UTILS_H

8
include/openmc/surface.h
View file

@ -32,13 +32,15 @@ extern "C" const int BC_PERIODIC;
// Global variables
//==============================================================================
extern "C" int32_t n_surfaces;
class Surface;
extern std::vector<Surface*> surfaces;
namespace model {
extern std::vector<Surface*> surfaces;
extern std::map<int, int> surface_map;
} // namespace model
//==============================================================================
//! Coordinates for an axis-aligned cube that bounds a geometric object.
//==============================================================================

11
include/openmc/tallies/filter.h
View file

@ -29,7 +29,7 @@ public:
// Without an explicit instantiation of vector<FilterMatch>, the Intel compiler
// will complain about the threadprivate directive on filter_matches. Note that
// this has to happen *outside* of the openmc namespace
template class std::vector<openmc::FilterMatch>;
extern template class std::vector<openmc::FilterMatch>;
namespace openmc {
@ -77,13 +77,20 @@ public:
// Global variables
//==============================================================================
extern "C" int32_t n_filters;
namespace simulation {
extern std::vector<FilterMatch> filter_matches;
#pragma omp threadprivate(filter_matches)
} // namespace simulation
namespace model {
extern "C" int32_t n_filters;
extern std::vector<std::unique_ptr<Filter>> tally_filters;
} // namespace model
//==============================================================================
extern "C" void free_memory_tally_c();

36
include/openmc/timer.h
View file

@ -5,6 +5,27 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
class Timer;
namespace simulation {
extern Timer time_active;
extern Timer time_bank;
extern Timer time_bank_sample;
extern Timer time_bank_sendrecv;
extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;
} // namespace simulation
//==============================================================================
//! Class for measuring time elapsed
//==============================================================================
@ -34,21 +55,6 @@ private:
double elapsed_ {0.0}; //!< elasped time in [s]
};
//==============================================================================
// Global variables
//==============================================================================
extern Timer time_active;
extern Timer time_bank;
extern Timer time_bank_sample;
extern Timer time_bank_sendrecv;
extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;
//==============================================================================
// Non-member functions
//==============================================================================

5
man/man1/openmc.1
View file

@ -53,11 +53,6 @@ to locate HDF5 format cross section libraries if the user has not specified the
Indicates the default path to an HDF5 file that contains multi-group cross
section libraries if the user has not specified the <cross_sections> tag in
.I materials.xml\fP.
.TP
.B OPENMC_MULTIPOLE_LIBRARY
Indicates the default path to a directory containing windowed multipole data if
the user has not specified the <multipole_library> tag in
.I materials.xml\fP.
.SH LICENSE
Copyright \(co 2011-2018 Massachusetts Institute of Technology and OpenMC
contributors.

14
openmc/_utils.py
View file

@ -1,3 +1,4 @@
import hashlib
import os.path
from pathlib import Path
from urllib.parse import urlparse
@ -6,13 +7,15 @@ from urllib.request import urlopen
_BLOCK_SIZE = 16384
def download(url):
def download(url, checksum=None):
"""Download file from a URL
Parameters
----------
url : str
URL from which to download
checksum : str or None
MD5 checksum to check against
Returns
-------
@ -46,4 +49,13 @@ def download(url):
downloaded, downloaded * 100. / file_size)
print(status + '\b'*len(status), end='')
print('')
if checksum is not None:
downloadsum = hashlib.md5(open(basename, 'rb').read()).hexdigest()
if downloadsum != checksum:
raise IOError("MD5 checksum for {} does not match. If this is your first "
"time receiving this message, please re-run the script. "
"Otherwise, please contact OpenMC developers by emailing "
"openmc-users@googlegroups.com.".format(basename))
return basename

2
openmc/capi/core.py
View file

@ -52,7 +52,7 @@ _dll.openmc_simulation_init.restype = c_int
_dll.openmc_simulation_init.errcheck = _error_handler
_dll.openmc_simulation_finalize.restype = c_int
_dll.openmc_simulation_finalize.errcheck = _error_handler
_dll.openmc_statepoint_write.argtypes = [POINTER(c_char_p), POINTER(c_bool)]
_dll.openmc_statepoint_write.argtypes = [c_char_p, POINTER(c_bool)]
_dll.openmc_statepoint_write.restype = c_int
_dll.openmc_statepoint_write.errcheck = _error_handler

10
openmc/capi/material.py
View file

@ -35,7 +35,7 @@ _dll.openmc_material_get_densities.errcheck = _error_handler
_dll.openmc_material_get_volume.argtypes = [c_int32, POINTER(c_double)]
_dll.openmc_material_get_volume.restype = c_int
_dll.openmc_material_get_volume.errcheck = _error_handler
_dll.openmc_material_set_density.argtypes = [c_int32, c_double]
_dll.openmc_material_set_density.argtypes = [c_int32, c_double, c_char_p]
_dll.openmc_material_set_density.restype = c_int
_dll.openmc_material_set_density.errcheck = _error_handler
_dll.openmc_material_set_densities.argtypes = [
@ -178,16 +178,18 @@ class Material(_FortranObjectWithID):
"""
_dll.openmc_material_add_nuclide(self._index, name.encode(), density)
def set_density(self, density):
def set_density(self, density, units='atom/b-cm'):
"""Set density of a material.
Parameters
----------
density : float
Density in atom/b-cm
Density
units : {'atom/b-cm', 'g/cm3'}
Units for density
"""
_dll.openmc_material_set_density(self._index, density)
_dll.openmc_material_set_density(self._index, density, units.encode())
def set_densities(self, nuclides, densities):
"""Set the densities of a list of nuclides in a material

2
openmc/data/__init__.py
View file

@ -1,5 +1,5 @@
# Version of HDF5 nuclear data format
HDF5_VERSION_MAJOR = 1
HDF5_VERSION_MAJOR = 2
HDF5_VERSION_MINOR = 0
HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)

7
openmc/data/ace.py
View file

@ -106,7 +106,7 @@ def ascii_to_binary(ascii_file, binary_file):
"""
# Open ASCII file
ascii = open(ascii_file, 'r')
ascii = open(str(ascii_file), 'r')
# Set default record length
record_length = 4096
@ -116,7 +116,7 @@ def ascii_to_binary(ascii_file, binary_file):
ascii.close()
# Open binary file
binary = open(binary_file, 'wb')
binary = open(str(binary_file), 'wb')
idx = 0
@ -228,8 +228,9 @@ class Library(EqualityMixin):
self.tables = []
# Determine whether file is ASCII or binary
filename = str(filename)
try:
fh = open(str(filename), 'rb')
fh = open(filename, 'rb')
# Grab 10 lines of the library
sb = b''.join([fh.readline() for i in range(10)])

2
openmc/data/endf.py
View file

@ -348,7 +348,7 @@ def get_evaluations(filename):
"""
evaluations = []
with open(filename, 'r') as fh:
with open(str(filename), 'r') as fh:
while True:
pos = fh.tell()
line = fh.readline()

4
openmc/data/fission_energy.py
View file

@ -279,11 +279,11 @@ class FissionEnergyRelease(EqualityMixin):
# the delayed neutron fraction is so small that the difference
# is negligible. MT=18 (n, fission) might not be available so
# try MT=19 (n, f) as well.
if 18 in incident_neutron.reactions:
if 18 in incident_neutron and not incident_neutron[18].redundant:
nu = [p.yield_ for p in incident_neutron[18].products
if p.particle == 'neutron'
and p.emission_mode in ('prompt', 'total')]
elif 19 in incident_neutron.reactions:
elif 19 in incident_neutron:
nu = [p.yield_ for p in incident_neutron[19].products
if p.particle == 'neutron'
and p.emission_mode in ('prompt', 'total')]

2
openmc/data/function.py
View file

@ -518,7 +518,7 @@ class Sum(EqualityMixin):
"""
def __init__(self, functions):
self.functions = functions
self.functions = list(functions)
def __call__(self, x):
return sum(f(x) for f in self.functions)

24
openmc/data/library.py
View file

@ -52,17 +52,13 @@ class DataLibrary(EqualityMixin):
Path to the file to be registered.
"""
with h5py.File(filename, 'r') as h5file:
# Support pathlib
# TODO: Remove when support is Python 3.6+ only
filename = str(filename)
materials = []
if 'filetype' in h5file.attrs:
filetype = h5file.attrs['filetype'].decode().lstrip('data_')
else:
filetype = 'neutron'
for name in h5file:
if name.startswith('c_'):
filetype = 'thermal'
materials.append(name)
with h5py.File(filename, 'r') as h5file:
filetype = h5file.attrs['filetype'].decode()[5:]
materials = list(h5file)
library = {'path': filename, 'type': filetype, 'materials': materials}
self.libraries.append(library)
@ -84,7 +80,7 @@ class DataLibrary(EqualityMixin):
if common_dir == '':
common_dir = '.'
if os.path.relpath(common_dir, os.path.dirname(path)) != '.':
if os.path.relpath(common_dir, os.path.dirname(str(path))) != '.':
dir_element = ET.SubElement(root, "directory")
dir_element.text = os.path.realpath(common_dir)
@ -99,7 +95,7 @@ class DataLibrary(EqualityMixin):
# Write XML file
tree = ET.ElementTree(root)
tree.write(path, xml_declaration=True, encoding='utf-8',
tree.write(str(path), xml_declaration=True, encoding='utf-8',
method='xml')
@classmethod
@ -131,7 +127,9 @@ class DataLibrary(EqualityMixin):
raise ValueError("Either path or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
check_type('path', path, str)
# Convert to string to support pathlib
# TODO: Remove when support is Python 3.6+ only
path = str(path)
tree = ET.parse(path)
root = tree.getroot()

4
openmc/data/multipole.py
View file

@ -333,7 +333,7 @@ class WindowedMultipole(EqualityMixin):
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(group_or_filename, 'r')
h5file = h5py.File(str(group_or_filename), 'r')
# Make sure version matches
if 'version' in h5file.attrs:
@ -515,7 +515,7 @@ class WindowedMultipole(EqualityMixin):
"""
# Open file and write version.
with h5py.File(path, mode, libver=libver) as f:
with h5py.File(str(path), mode, libver=libver) as f:
f.attrs['filetype'] = np.string_('data_wmp')
f.attrs['version'] = np.array(WMP_VERSION)

150
openmc/data/neutron.py
View file

@ -2,7 +2,6 @@ import sys
from collections import OrderedDict
from collections.abc import Iterable, Mapping, MutableMapping
from io import StringIO
from itertools import chain
from math import log10
from numbers import Integral, Real
import os
@ -86,9 +85,6 @@ class IncidentNeutron(EqualityMixin):
Resonance parameters
resonance_covariance : openmc.data.ResonanceCovariance or None
Covariance for resonance parameters
redundant_reactions : collections.OrderedDict
Contains redundant cross sections, e.g., the total cross section. The keys
are the MT values and the values are Reaction objects.
temperatures : list of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings of the temperature,
@ -113,18 +109,15 @@ class IncidentNeutron(EqualityMixin):
self.energy = {}
self._fission_energy = None
self.reactions = OrderedDict()
self.redundant_reactions = OrderedDict()
self._urr = {}
self._resonances = None
def __contains__(self, mt):
return mt in self.reactions or mt in self.redundant_reactions
return mt in self.reactions
def __getitem__(self, mt):
if mt in self.reactions:
return self.reactions[mt]
elif mt in self.redundant_reactions:
return self.redundant_reactions[mt]
else:
raise KeyError('No reaction with MT={}.'.format(mt))
@ -170,10 +163,6 @@ class IncidentNeutron(EqualityMixin):
def resonance_covariance(self):
return self._resonance_covariance
@property
def redundant_reactions(self):
return self._redundant_reactions
@property
def urr(self):
return self._urr
@ -237,11 +226,6 @@ class IncidentNeutron(EqualityMixin):
res_cov.ResonanceCovariances)
self._resonance_covariance = resonance_covariance
@redundant_reactions.setter
def redundant_reactions(self, redundant_reactions):
cv.check_type('redundant reactions', redundant_reactions, Mapping)
self._redundant_reactions = redundant_reactions
@urr.setter
def urr(self, urr):
cv.check_type('probability table dictionary', urr, MutableMapping)
@ -288,7 +272,7 @@ class IncidentNeutron(EqualityMixin):
self.energy[strT] = data.energy[strT]
# Add normal and redundant reactions
for mt in chain(data.reactions, data.redundant_reactions):
for mt in data.reactions:
if mt in self:
self[mt].xs[strT] = data[mt].xs[strT]
else:
@ -406,24 +390,14 @@ class IncidentNeutron(EqualityMixin):
have cross sections provided.
"""
if mt in self.reactions:
return [mt]
elif mt in SUM_RULES:
mts = SUM_RULES[mt]
mts = []
if mt in SUM_RULES:
for mt_i in SUM_RULES[mt]:
mts += self.get_reaction_components(mt_i)
if mts:
return mts
else:
return []
complete = False
while not complete:
new_mts = []
complete = True
for i, mt_i in enumerate(mts):
if mt_i in self.reactions:
new_mts.append(mt_i)
elif mt_i in SUM_RULES:
new_mts += SUM_RULES[mt_i]
complete = False
mts = new_mts
return mts
return [mt] if mt in self else []
def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export incident neutron data to an HDF5 file.
@ -446,7 +420,7 @@ class IncidentNeutron(EqualityMixin):
'originated from an ENDF file.')
# Open file and write version
f = h5py.File(path, mode, libver=libver)
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_neutron')
f.attrs['version'] = np.array(HDF5_VERSION)
@ -472,6 +446,15 @@ class IncidentNeutron(EqualityMixin):
# Write reaction data
rxs_group = g.create_group('reactions')
for rx in self.reactions.values():
# Skip writing redundant reaction if it doesn't have photon
# production or is a summed transmutation reaction. MT=4 is also
# sometimes needed for probability tables.
if rx.redundant:
photon_rx = any(p.particle == 'photon' for p in rx.products)
transmutation_rx = (rx.mt in (16, 103, 104, 105, 106, 107))
if not (photon_rx or transmutation_rx or rx.mt == 4):
continue
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
rx.to_hdf5(rx_group)
@ -480,12 +463,6 @@ class IncidentNeutron(EqualityMixin):
tgroup = g.create_group('total_nu')
rx.derived_products[0].to_hdf5(tgroup)
# Write redundant reaction data only for reactions with photon production
for rx in self.redundant_reactions.values():
if any(p.particle == 'photon' for p in rx.products):
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
rx.to_hdf5(rx_group)
# Write unresolved resonance probability tables
if self.urr:
urr_group = g.create_group('urr')
@ -520,16 +497,12 @@ class IncidentNeutron(EqualityMixin):
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(group_or_filename, 'r')
h5file = h5py.File(str(group_or_filename), 'r')
# Make sure version matches
if 'version' in h5file.attrs:
major, minor = h5file.attrs['version']
if major != HDF5_VERSION_MAJOR:
raise IOError(
'HDF5 data format uses version {}.{} whereas your '
'installation of the OpenMC Python API expects version '
'{}.x.'.format(major, minor, HDF5_VERSION_MAJOR))
# For now all versions of HDF5 data can be read
else:
raise IOError(
'HDF5 data does not indicate a version. Your installation of '
@ -561,10 +534,7 @@ class IncidentNeutron(EqualityMixin):
for name, obj in sorted(rxs_group.items()):
if name.startswith('reaction_'):
rx = Reaction.from_hdf5(obj, data.energy)
if rx.redundant:
data.redundant_reactions[rx.mt] = rx
else:
data.reactions[rx.mt] = rx
data.reactions[rx.mt] = rx
# Read total nu data if available
if rx.mt in (18, 19, 20, 21, 38) and 'total_nu' in group:
@ -577,7 +547,8 @@ class IncidentNeutron(EqualityMixin):
if mt_sum not in data:
rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data]
if len(rxs) > 0:
data.redundant_reactions[mt_sum] = rx = Reaction(mt_sum)
data.reactions[mt_sum] = rx = Reaction(mt_sum)
rx.redundant = True
if rx.mt == 18 and 'total_nu' in group:
tgroup = group['total_nu']
rx.derived_products.append(Product.from_hdf5(tgroup))
@ -654,18 +625,21 @@ class IncidentNeutron(EqualityMixin):
total = Reaction(1)
total.xs[strT] = Tabulated1D(energy, total_xs)
total.redundant = True
data.redundant_reactions[1] = total
data.reactions[1] = total
if np.count_nonzero(absorption_xs) > 0:
absorption = Reaction(27)
absorption = Reaction(101)
absorption.xs[strT] = Tabulated1D(energy, absorption_xs)
absorption.redundant = True
data.redundant_reactions[27] = absorption
data.reactions[101] = absorption
# Read each reaction
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
rx = Reaction.from_ace(ace, i)
# Don't include gas production / damage cross sections
if 200 < rx.mt < 219 or rx.mt == 444:
continue
data.reactions[rx.mt] = rx
# Some photon production reactions may be assigned to MTs that don't
@ -683,23 +657,39 @@ class IncidentNeutron(EqualityMixin):
continue
# Create redundant reaction with appropriate cross section
rx = Reaction(mt)
mts = data.get_reaction_components(mt)
if len(mts) == 0:
warn('Photon production is present for MT={} but no '
'reaction components exist.'.format(mt))
continue
xss = [data.reactions[mt_i].xs[strT] for mt_i in mts]
idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
else 0 for xs in xss])
rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
rx.xs[strT]._threshold_idx = idx
rx.redundant = True
# Determine redundant cross section
rx = data._get_redundant_reaction(mt, mts)
rx.products += _get_photon_products_ace(ace, rx)
data.reactions[mt] = rx
# For transmutation reactions, sometimes only individual levels are
# present in an ACE file, e.g. MT=600-649 instead of the summation
# MT=103. In this case, if a user wants to tally (n,p), OpenMC doesn't
# know about the total cross section. Here, we explicitly create a
# redundant reaction for this purpose.
for mt in (16, 103, 104, 105, 106, 107):
if mt not in data:
# Determine if any individual levels are present
mts = data.get_reaction_components(mt)
if len(mts) == 0:
continue
# Determine redundant cross section
rx.products += _get_photon_products_ace(ace, rx)
data.redundant_reactions[mt] = rx
rx = data._get_redundant_reaction(mt, mts)
data.reactions[mt] = rx
# Make sure redundant cross sections that are present in an ACE file get
# marked as such
for rx in data:
mts = data.get_reaction_components(rx.mt)
if mts != [rx.mt]:
rx.redundant = True
# Read unresolved resonance probability tables
urr = ProbabilityTables.from_ace(ace)
@ -842,3 +832,33 @@ class IncidentNeutron(EqualityMixin):
data[2].xs['0K'] = xs
return data
def _get_redundant_reaction(self, mt, mts):
"""Create redundant reaction from its components
Parameters
----------
mt : int
MT value of the desired reaction
mts : iterable of int
MT values of its components
Returns
-------
openmc.Reaction
Redundant reaction
"""
# Get energy grid
strT = self.temperatures[0]
energy = self.energy[strT]
rx = Reaction(mt)
xss = [self.reactions[mt_i].xs[strT] for mt_i in mts]
idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
else 0 for xs in xss])
rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
rx.xs[strT]._threshold_idx = idx
rx.redundant = True
return rx

8
openmc/data/njoy.py
View file

@ -146,14 +146,14 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
"""
if input_filename is not None:
with open(input_filename, 'w') as f:
with open(str(input_filename), 'w') as f:
f.write(commands)
with tempfile.TemporaryDirectory() as tmpdir:
# Copy evaluations to appropriates 'tapes'
for tape_num, filename in tapein.items():
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
shutil.copy(filename, tmpfilename)
shutil.copy(str(filename), tmpfilename)
# Start up NJOY process
njoy = Popen([njoy_exec], cwd=tmpdir, stdin=PIPE, stdout=PIPE,
@ -182,7 +182,7 @@ def run(commands, tapein, tapeout, input_filename=None, stdout=False,
for tape_num, filename in tapeout.items():
tmpfilename = os.path.join(tmpdir, 'tape{}'.format(tape_num))
if os.path.isfile(tmpfilename):
shutil.move(tmpfilename, filename)
shutil.move(tmpfilename, str(filename))
def make_pendf(filename, pendf='pendf', error=0.001, stdout=False):
@ -422,7 +422,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
commands = ""
nendf, nthermal_endf, npendf = 20, 21, 22
tapein = {nendf: filename, nthermal_endf:filename_thermal}
tapein = {nendf: filename, nthermal_endf: filename_thermal}
tapeout = {}
# reconr

10
openmc/data/photon.py
View file

@ -373,9 +373,6 @@ class IncidentPhoton(EqualityMixin):
excitation energy), 's_collision' (collision stopping power in
[eV cm\ :sup:`2`/g]), and 's_radiative' (radiative stopping power in
[eV cm\ :sup:`2`/g])
redundant_reactions : collections.OrderedDict
Contains redundant cross sections. The keys are MT values and the values
are instances of :class:`PhotonReaction`.
"""
@ -383,19 +380,16 @@ class IncidentPhoton(EqualityMixin):
self.atomic_number = atomic_number
self._atomic_relaxation = None
self.reactions = OrderedDict()
self.redundant_reactions = OrderedDict()
self.compton_profiles = {}
self.stopping_powers = {}
self.bremsstrahlung = {}
def __contains__(self, mt):
return mt in self.reactions or mt in self.redundant_reactions
return mt in self.reactions
def __getitem__(self, mt):
if mt in self.reactions:
return self.reactions[mt]
elif mt in self.redundant_reactions:
return self.redundant_reactions[mt]
else:
raise KeyError('No reaction with MT={}.'.format(mt))
@ -668,7 +662,7 @@ class IncidentPhoton(EqualityMixin):
"""
# Open file and write version
f = h5py.File(path, mode, libver=libver)
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_photon')
if 'version' not in f.attrs:
f.attrs['version'] = np.array(HDF5_VERSION)

4
openmc/data/thermal.py
View file

@ -278,7 +278,7 @@ class ThermalScattering(EqualityMixin):
"""
# Open file and write version
f = h5py.File(path, mode, libver=libver)
f = h5py.File(str(path), mode, libver=libver)
f.attrs['filetype'] = np.string_('data_thermal')
f.attrs['version'] = np.array(HDF5_VERSION)
@ -387,7 +387,7 @@ class ThermalScattering(EqualityMixin):
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(group_or_filename, 'r')
h5file = h5py.File(str(group_or_filename), 'r')
# Make sure version matches
if 'version' in h5file.attrs:

24
openmc/material.py
View file

@ -989,18 +989,12 @@ class Materials(cv.CheckedList):
continuous-energy calculations and
:envvar:`OPENMC_MG_CROSS_SECTIONS` will be used for multi-group
calculations to find the path to the HDF5 cross section file.
multipole_library : str
Indicates the path to a directory containing a windowed multipole
cross section library. If it is not set, the
:envvar:`OPENMC_MULTIPOLE_LIBRARY` environment variable will be used. A
multipole library is optional.
"""
def __init__(self, materials=None):
super().__init__(Material, 'materials collection')
self._cross_sections = None
self._multipole_library = None
if materials is not None:
self += materials
@ -1009,20 +1003,11 @@ class Materials(cv.CheckedList):
def cross_sections(self):
return self._cross_sections
@property
def multipole_library(self):
return self._multipole_library
@cross_sections.setter
def cross_sections(self, cross_sections):
cv.check_type('cross sections', cross_sections, str)
self._cross_sections = cross_sections
@multipole_library.setter
def multipole_library(self, multipole_library):
cv.check_type('cross sections', multipole_library, str)
self._multipole_library = multipole_library
def append(self, material):
"""Append material to collection
@ -1060,11 +1045,6 @@ class Materials(cv.CheckedList):
element = ET.SubElement(root_element, "cross_sections")
element.text = str(self._cross_sections)
def _create_multipole_library_subelement(self, root_element):
if self._multipole_library is not None:
element = ET.SubElement(root_element, "multipole_library")
element.text = str(self._multipole_library)
def export_to_xml(self, path='materials.xml'):
"""Export material collection to an XML file.
@ -1077,7 +1057,6 @@ class Materials(cv.CheckedList):
root_element = ET.Element("materials")
self._create_cross_sections_subelement(root_element)
self._create_multipole_library_subelement(root_element)
self._create_material_subelements(root_element)
# Clean the indentation in the file to be user-readable
@ -1114,8 +1093,5 @@ class Materials(cv.CheckedList):
xs = tree.find('cross_sections')
if xs is not None:
materials.cross_sections = xs.text
mpl = tree.find('multipole_library')
if mpl is not None:
materials.multipole_library = mpl.text
return materials

891
openmc/model/triso.py
View file

File diff suppressed because it is too large Load diff

2
openmc/region.py
View file

@ -475,7 +475,7 @@ class Complement(Region):
>>> xr = openmc.XPlane(x0=10.0)
>>> yl = openmc.YPlane(y0=-10.0)
>>> yr = openmc.YPlane(y0=10.0)
>>> inside_box = +xl & -xr & +yl & -yl
>>> inside_box = +xl & -xr & +yl & -yr
>>> outside_box = ~inside_box
>>> type(outside_box)
<class 'openmc.region.Complement'>

1
openmc/settings.py
View file

@ -180,7 +180,6 @@ class Settings(object):
self._confidence_intervals = None
self._cross_sections = None
self._electron_treatment = None
self._multipole_library = None
self._photon_transport = None
self._ptables = None
self._run_cmfd = None

4
openmc/surface.py
View file

@ -62,7 +62,7 @@ class Surface(IDManagerMixin):
self.boundary_type = boundary_type
# A dictionary of the quadratic surface coefficients
# Key - coefficeint name
# Key - coefficient name
# Value - coefficient value
self._coefficients = {}
@ -1661,7 +1661,7 @@ class Quadric(Surface):
a, b, c, d, e, f, g, h, j, k : float, optional
coefficients for the surface. All default to 0.
name : str, optional
Name of the sphere. If not specified, the name will be the empty string.
Name of the surface. If not specified, the name will be the empty string.
Attributes
----------

99
scripts/openmc-convert-lib80x-data Executable file
View file

@ -0,0 +1,99 @@
#!/usr/bin/env python3
import argparse
from collections import defaultdict
import glob
import os
import openmc.data
description = """
Convert ENDF/B-VIII.0 ACE data from LANL into an HDF5 library
that can be used by OpenMC. This assumes that you have a directory containing
subdirectories 'Lib80x' and 'ENDF80SaB'.
"""
class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
argparse.RawDescriptionHelpFormatter):
pass
parser = argparse.ArgumentParser(
description=description,
formatter_class=CustomFormatter
)
parser.add_argument('-o', '--output_dir', default='lib80x_hdf5',
help='Directory to create new library in')
parser.add_argument('--libver', choices=['earliest', 'latest'],
default='earliest', help="Output HDF5 versioning. Use "
"'earliest' for backwards compatibility or 'latest' for "
"performance")
parser.add_argument('--datadir', help='Directory containing Lib80x and ENDF80SaB',
default=os.curdir)
args = parser.parse_args()
assert os.path.isdir(args.datadir)
# Get a list of all ACE files
lib80x = glob.glob(os.path.join(args.datadir, 'Lib80x', '**', '*.80?nc'), recursive=True)
lib80sab = glob.glob(os.path.join(args.datadir, 'ENDF80SaB', '**', '*.??t'), recursive=True)
# Find and fix B10 ACE files
b10files = glob.glob(os.path.join(args.datadir, 'Lib80x', '**', '5010.80?nc'), recursive=True)
nxs1_position = 523
for filename in b10files:
with open(filename, 'r+') as fh:
# Read NXS(1)
fh.seek(nxs1_position)
nxs1 = int(fh.read(5))
# Increase length to match actual length of XSS, but make sure this
# isn't done twice by checking the current length
if nxs1 < 86870:
fh.seek(nxs1_position)
fh.write(str(nxs1 + 53))
# Group together tables for the same nuclide
suffixes = defaultdict(list)
for filename in sorted(lib80x + lib80sab):
dirname, basename = os.path.split(filename)
zaid, xs = basename.split('.')
suffixes[os.path.join(dirname, zaid)].append(xs)
# Create output directory if it doesn't exist
if not os.path.isdir(args.output_dir):
os.mkdir(args.output_dir)
library = openmc.data.DataLibrary()
for basename, xs_list in sorted(suffixes.items()):
# Convert first temperature for the table
filename = '.'.join((basename, xs_list[0]))
print('Converting: ' + filename)
if filename.endswith('t'):
data = openmc.data.ThermalScattering.from_ace(filename)
else:
data = openmc.data.IncidentNeutron.from_ace(filename, 'mcnp')
# For each higher temperature, add cross sections to the existing table
for xs in xs_list[1:]:
filename = '.'.join((basename, xs))
print('Adding: ' + filename)
if filename.endswith('t'):
data.add_temperature_from_ace(filename)
else:
data.add_temperature_from_ace(filename, 'mcnp')
# Export HDF5 file
h5_file = os.path.join(args.output_dir, data.name + '.h5')
print('Writing {}...'.format(h5_file))
data.export_to_hdf5(h5_file, 'w', libver=args.libver)
# Register with library
library.register_file(h5_file)
# Write cross_sections.xml
libpath = os.path.join(args.output_dir, 'cross_sections.xml')
library.export_to_xml(libpath)

120
scripts/openmc-get-multipole-data
View file

@ -1,120 +0,0 @@
#!/usr/bin/env python3
import os
import shutil
import subprocess
import sys
import tarfile
import glob
import hashlib
import argparse
from urllib.request import urlopen
description = """
Download and extract windowed multipole data based on ENDF/B-VII.1.
"""
class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter,
argparse.RawDescriptionHelpFormatter):
pass
parser = argparse.ArgumentParser(
description=description,
formatter_class=CustomFormatter
)
parser.add_argument('-b', '--batch', action='store_true',
help='supresses standard in')
args = parser.parse_args()
baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/'
files = ['WMP_Library_v1.1.tar.gz']
checksums = ['8523895928dd6ba63fba803e3a45d4f3']
block_size = 16384
# ==============================================================================
# DOWNLOAD FILES FROM GITHUB REPO
filesComplete = []
for f in files:
# Establish connection to URL
url = baseUrl + f
req = urlopen(url)
# Get file size from header
if sys.version_info[0] < 3:
file_size = int(req.info().getheaders('Content-Length')[0])
else:
file_size = req.length
downloaded = 0
# Remove GitHub junk from the file name.
fname = f[:-9] if f.endswith('?raw=true') else f
# Check if file already downloaded
if os.path.exists(fname):
if os.path.getsize(fname) == file_size:
print('Skipping ' + fname)
filesComplete.append(fname)
continue
else:
overwrite = input('Overwrite {0}? ([y]/n) '.format(fname))
if overwrite.lower().startswith('n'):
continue
# Copy file to disk
print('Downloading {0}... '.format(f), end='')
with open(fname, 'wb') as fh:
while True:
chunk = req.read(block_size)
if not chunk: break
fh.write(chunk)
downloaded += len(chunk)
status = '{0:10} [{1:3.2f}%]'.format(downloaded, downloaded * 100. / file_size)
print(status + chr(8)*len(status), end='')
print('')
filesComplete.append(fname)
# ==============================================================================
# VERIFY MD5 CHECKSUMS
print('Verifying MD5 checksums...')
for f, checksum in zip(files, checksums):
fname = f[:-9] if f.endswith('?raw=true') else f
downloadsum = hashlib.md5(open(fname, 'rb').read()).hexdigest()
if downloadsum != checksum:
raise IOError("MD5 checksum for {} does not match. If this is your first "
"time receiving this message, please re-run the script. "
"Otherwise, please contact OpenMC developers by emailing "
"openmc-users@googlegroups.com.".format(f))
# ==============================================================================
# EXTRACT FILES FROM TGZ
for f in files:
fname = f[:-9] if f.endswith('?raw=true') else f
if fname not in filesComplete:
continue
# Extract files
with tarfile.open(fname, 'r') as tgz:
print('Extracting {0}...'.format(fname))
tgz.extractall(path='')
# ==============================================================================
# PROMPT USER TO DELETE .TAR.GZ FILES
# Ask user to delete
if not args.batch:
response = input('Delete *.tar.gz files? ([y]/n) ')
else:
response = 'y'
# Delete files if requested
if not response or response.lower().startswith('y'):
for f in files:
if os.path.exists(f):
print('Removing {0}...'.format(f))
os.remove(f)

164
scripts/openmc-make-test-data Executable file
View file

@ -0,0 +1,164 @@
#!/usr/bin/env python3
"""
Download ENDF/B-VII.1 ENDF and ACE files from NNDC and WMP files from GitHub and
generate a full HDF5 library with incident neutron, incident photon, thermal
scattering data, and windowed multipole data. This data is used for OpenMC's
regression test suite.
"""
import glob
import os
from pathlib import Path
import tarfile
import tempfile
from urllib.parse import urljoin
import zipfile
import openmc.data
from openmc._utils import download
base_ace = 'http://www.nndc.bnl.gov/endf/b7.1/aceFiles/'
base_endf = 'http://www.nndc.bnl.gov/endf/b7.1/zips/'
base_wmp = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/'
files = [
(base_ace, 'ENDF-B-VII.1-neutron-293.6K.tar.gz', '9729a17eb62b75f285d8a7628ace1449'),
(base_ace, 'ENDF-B-VII.1-tsl.tar.gz', 'e17d827c92940a30f22f096d910ea186'),
(base_endf, 'ENDF-B-VII.1-neutrons.zip', 'e5d7f441fc4c92893322c24d1725e29c'),
(base_endf, 'ENDF-B-VII.1-photoat.zip', '5192f94e61f0b385cf536f448ffab4a4'),
(base_endf, 'ENDF-B-VII.1-atomic_relax.zip', 'fddb6035e7f2b6931e51a58fc754bd10'),
(base_wmp, 'WMP_Library_v1.1.tar.gz', '8523895928dd6ba63fba803e3a45d4f3')
]
def fix_zaid(table, old, new):
filename = os.path.join('tsl', table)
with open(filename, 'r') as fh:
text = fh.read()
text = text.replace(old, new, 1)
with open(filename, 'w') as fh:
fh.write(text)
pwd = Path.cwd()
output_dir = pwd / 'nndc_hdf5'
os.makedirs('nndc_hdf5/photon', exist_ok=True)
with tempfile.TemporaryDirectory() as tmpdir:
# Temporarily change dir
os.chdir(tmpdir)
# =========================================================================
# Download files from NNDC server
for base, fname, checksum in files:
download(urljoin(base, fname), checksum)
# =========================================================================
# EXTRACT FILES FROM TGZ
for _, f, _ in files:
print('Extracting {}...'.format(f))
path = Path(f)
if path.suffix == '.gz':
with tarfile.open(f, 'r') as tgz:
if 'tsl' in f:
tgz.extractall(path='tsl')
else:
tgz.extractall()
elif path.suffix == '.zip':
zipfile.ZipFile(f).extractall()
# =========================================================================
# FIX ZAID ASSIGNMENTS FOR VARIOUS S(A,B) TABLES
print('Fixing ZAIDs for S(a,b) tables')
fix_zaid('bebeo.acer', '8016', ' 0')
fix_zaid('obeo.acer', '4009', ' 0')
library = openmc.data.DataLibrary()
# =========================================================================
# INCIDENT NEUTRON DATA
neutron_files = sorted(glob.glob('ENDF-B-VII.1-neutron-293.6K/*.ace'))
for f in neutron_files:
print('Converting {}...'.format(os.path.basename(f)))
data = openmc.data.IncidentNeutron.from_ace(f)
# Check for fission energy release data
endf_filename = 'neutrons/n-{:03}_{}_{:03}{}.endf'.format(
data.atomic_number,
data.atomic_symbol,
data.mass_number,
'm{}'.format(data.metastable) if data.metastable else ''
)
ev = openmc.data.endf.Evaluation(endf_filename)
if (1, 458) in ev.section:
endf_data = openmc.data.IncidentNeutron.from_endf(ev)
data.fission_energy = endf_data.fission_energy
# Add 0K elastic scattering data for select nuclides
if data.name in ('U235', 'U238', 'Pu239'):
data.add_elastic_0K_from_endf(endf_filename)
# Determine filename
outfile = output_dir / (data.name + '.h5')
data.export_to_hdf5(outfile, 'w', 'earliest')
# Register with library
library.register_file(outfile)
# =========================================================================
# THERMAL SCATTERING DATA
thermal_files = sorted(glob.glob('tsl/*.acer'))
for f in thermal_files:
print('Converting {}...'.format(os.path.basename(f)))
data = openmc.data.ThermalScattering.from_ace(f)
# Determine filename
outfile = output_dir / (data.name + '.h5')
data.export_to_hdf5(outfile, 'w', 'earliest')
# Register with library
library.register_file(outfile)
# =========================================================================
# INCIDENT PHOTON DATA
for z in range(1, 101):
element = openmc.data.ATOMIC_SYMBOL[z]
print('Generating HDF5 file for Z={} ({})...'.format(z, element))
# Generate instance of IncidentPhoton
photo_file = Path('photoat') / 'photoat-{:03}_{}_000.endf'.format(z, element)
atom_file = Path('atomic_relax') / 'atom-{:03}_{}_000.endf'.format(z, element)
data = openmc.data.IncidentPhoton.from_endf(photo_file, atom_file)
# Write HDF5 file and register it
outfile = output_dir / 'photon' / (element + '.h5')
data.export_to_hdf5(outfile, 'w', 'earliest')
library.register_file(outfile)
# =========================================================================
# WINDOWED MULTIPOLE DATA
# Move data into output directory
os.rename('WMP_Library', str(output_dir / 'wmp'))
# Add multipole data to library
for f in sorted(glob.glob('{}/wmp/*.h5'.format(output_dir))):
print('Registering WMP file {}...'.format(f))
library.register_file(f)
library.export_to_xml(output_dir / 'cross_sections.xml')
# =========================================================================
# CREATE TARBALL AND MOVE BACK
print('Creating compressed archive...')
test_tar = pwd / 'nndc_hdf5_test.tar.xz'
with tarfile.open(str(test_tar), 'w:xz') as txz:
txz.add('nndc_hdf5')
# Change back to original directory
os.chdir(str(pwd))

3
src/api.F90
View file

@ -79,6 +79,9 @@ contains
subroutine free_memory_mesh() bind(C)
end subroutine free_memory_mesh
subroutine free_memory_settings() bind(C)
end subroutine free_memory_settings
end interface
call free_memory_geometry()

108
src/cell.cpp
View file

@ -21,6 +21,8 @@ namespace openmc {
// Global variables
//==============================================================================
namespace model {
int32_t n_cells {0};
std::vector<Cell*> cells;
@ -29,6 +31,8 @@ std::unordered_map<int32_t, int32_t> cell_map;
std::vector<Universe*> universes;
std::unordered_map<int32_t, int32_t> universe_map;
} // namespace model
//==============================================================================
//! Convert region specification string to integer tokens.
//!
@ -197,7 +201,7 @@ Universe::to_hdf5(hid_t universes_group) const
// Write the contained cells.
if (cells_.size() > 0) {
std::vector<int32_t> cell_ids;
for (auto i_cell : cells_) cell_ids.push_back(cells[i_cell]->id_);
for (auto i_cell : cells_) cell_ids.push_back(model::cells[i_cell]->id_);
write_dataset(group, "cells", cell_ids);
}
@ -209,7 +213,7 @@ Universe::to_hdf5(hid_t universes_group) const
//==============================================================================
CSGCell::CSGCell() {} // empty constructor
CSGCell::CSGCell(pugi::xml_node cell_node)
{
if (check_for_node(cell_node, "id")) {
@ -314,7 +318,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
// Convert user IDs to surface indices.
for (auto& r : region_) {
if (r < OP_UNION) {
r = copysign(surface_map[abs(r)] + 1, r);
r = copysign(model::surface_map[abs(r)] + 1, r);
}
}
@ -418,7 +422,7 @@ CSGCell::distance(Position r, Direction u, int32_t on_surface) const
// Calculate the distance to this surface.
// Note the off-by-one indexing
bool coincident {token == on_surface};
double d {surfaces[abs(token)-1]->distance(r, u, coincident)};
double d {model::surfaces[abs(token)-1]->distance(r, u, coincident)};
// Check if this distance is the new minimum.
if (d < min_dist) {
@ -446,7 +450,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
write_string(group, "name", name_, false);
}
write_dataset(group, "universe", universes[universe_]->id_);
write_dataset(group, "universe", model::universes[universe_]->id_);
// Write the region specification.
if (!region_.empty()) {
@ -464,7 +468,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
} else {
// Note the off-by-one indexing
region_spec << " "
<< copysign(surfaces[abs(token)-1]->id_, token);
<< copysign(model::surfaces[abs(token)-1]->id_, token);
}
}
write_string(group, "region", region_spec.str(), false);
@ -476,7 +480,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
std::vector<int32_t> mat_ids;
for (auto i_mat : material_) {
if (i_mat != MATERIAL_VOID) {
mat_ids.push_back(materials[i_mat]->id_);
mat_ids.push_back(model::materials[i_mat]->id_);
} else {
mat_ids.push_back(MATERIAL_VOID);
}
@ -494,7 +498,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
} else if (type_ == FILL_UNIVERSE) {
write_dataset(group, "fill_type", "universe");
write_dataset(group, "fill", universes[fill_]->id_);
write_dataset(group, "fill", model::universes[fill_]->id_);
if (translation_ != Position(0, 0, 0)) {
write_dataset(group, "translation", translation_);
}
@ -505,7 +509,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
} else if (type_ == FILL_LATTICE) {
write_dataset(group, "fill_type", "lattice");
write_dataset(group, "lattice", lattices[fill_]->id_);
write_dataset(group, "lattice", model::lattices[fill_]->id_);
}
close_group(group);
@ -527,7 +531,7 @@ CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const
return false;
} else {
// Note the off-by-one indexing
bool sense = surfaces[abs(token)-1]->sense(r, u);
bool sense = model::surfaces[abs(token)-1]->sense(r, u);
if (sense != (token > 0)) {return false;}
}
}
@ -569,7 +573,7 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const
stack[i_stack] = false;
} else {
// Note the off-by-one indexing
bool sense = surfaces[abs(token)-1]->sense(r, u);
bool sense = model::surfaces[abs(token)-1]->sense(r, u);
stack[i_stack] = (sense == (token > 0));
}
}
@ -609,10 +613,10 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const
} else { // indicate that particle is lost
surf_idx = -1;
}
return {dist, surf_idx};
}
bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
{
moab::ErrorCode rval;
@ -622,7 +626,7 @@ bool DAGCell::contains(Position r, Direction u, int32_t on_surface) const
double pnt[3] = {r.x, r.y, r.z};
double dir[3] = {u.x, u.y, u.z};
rval = dagmc_ptr_->point_in_volume(vol, pnt, result, dir);
MB_CHK_ERR_CONT(rval);
MB_CHK_ERR_CONT(rval);
return result;
}
@ -638,23 +642,23 @@ extern "C" void
read_cells(pugi::xml_node* node)
{
// Count the number of cells.
for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;}
if (n_cells == 0) {
for (pugi::xml_node cell_node: node->children("cell")) {model::n_cells++;}
if (model::n_cells == 0) {
fatal_error("No cells found in geometry.xml!");
}
// Loop over XML cell elements and populate the array.
cells.reserve(n_cells);
for (pugi::xml_node cell_node: node->children("cell")) {
cells.push_back(new CSGCell(cell_node));
model::cells.reserve(model::n_cells);
for (pugi::xml_node cell_node : node->children("cell")) {
model::cells.push_back(new CSGCell(cell_node));
}
// Fill the cell map.
for (int i = 0; i < cells.size(); i++) {
int32_t id = cells[i]->id_;
auto search = cell_map.find(id);
if (search == cell_map.end()) {
cell_map[id] = i;
for (int i = 0; i < model::cells.size(); i++) {
int32_t id = model::cells[i]->id_;
auto search = model::cell_map.find(id);
if (search == model::cell_map.end()) {
model::cell_map[id] = i;
} else {
std::stringstream err_msg;
err_msg << "Two or more cells use the same unique ID: " << id;
@ -663,22 +667,24 @@ read_cells(pugi::xml_node* node)
}
// Populate the Universe vector and map.
for (int i = 0; i < cells.size(); i++) {
int32_t uid = cells[i]->universe_;
auto it = universe_map.find(uid);
if (it == universe_map.end()) {
universes.push_back(new Universe());
universes.back()->id_ = uid;
universes.back()->cells_.push_back(i);
universe_map[uid] = universes.size() - 1;
for (int i = 0; i < model::cells.size(); i++) {
int32_t uid = model::cells[i]->universe_;
auto it = model::universe_map.find(uid);
if (it == model::universe_map.end()) {
model::universes.push_back(new Universe());
model::universes.back()->id_ = uid;
model::universes.back()->cells_.push_back(i);
model::universe_map[uid] = model::universes.size() - 1;
} else {
universes[it->second]->cells_.push_back(i);
model::universes[it->second]->cells_.push_back(i);
}
}
universes.shrink_to_fit();
model::universes.shrink_to_fit();
// Allocate the cell overlap count if necessary.
if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0);
if (settings::check_overlaps) {
model::overlap_check_count.resize(model::cells.size(), 0);
}
}
//==============================================================================
@ -688,9 +694,9 @@ read_cells(pugi::xml_node* node)
extern "C" int
openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
{
if (index >= 1 && index <= cells.size()) {
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*cells[index - 1]};
Cell& c {*model::cells[index - 1]};
*type = c.type_;
if (c.type_ == FILL_MATERIAL) {
*indices = c.material_.data();
@ -710,9 +716,9 @@ extern "C" int
openmc_cell_set_fill(int32_t index, int type, int32_t n,
const int32_t* indices)
{
if (index >= 1 && index <= cells.size()) {
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*cells[index - 1]};
Cell& c {*model::cells[index - 1]};
if (type == FILL_MATERIAL) {
c.type_ = FILL_MATERIAL;
c.material_.clear();
@ -720,7 +726,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
int i_mat = indices[i];
if (i_mat == MATERIAL_VOID) {
c.material_.push_back(MATERIAL_VOID);
} else if (i_mat >= 1 && i_mat <= materials.size()) {
} else if (i_mat >= 1 && i_mat <= model::materials.size()) {
//TODO: off-by-one
c.material_.push_back(i_mat - 1);
} else {
@ -745,9 +751,9 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
extern "C" int
openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
{
if (index >= 1 && index <= cells.size()) {
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*cells[index - 1]};
Cell& c {*model::cells[index - 1]};
if (instance) {
if (*instance >= 0 && *instance < c.sqrtkT_.size()) {
@ -775,7 +781,7 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
//==============================================================================
extern "C" {
Cell* cell_pointer(int32_t cell_ind) {return cells[cell_ind];}
Cell* cell_pointer(int32_t cell_ind) {return model::cells[cell_ind];}
int32_t cell_id(Cell* c) {return c->id_;}
@ -785,9 +791,9 @@ extern "C" {
c->id_ = id;
// Find the index of this cell and update the cell map.
for (int i = 0; i < cells.size(); i++) {
if (cells[i] == c) {
cell_map[id] = i;
for (int i = 0; i < model::cells.size(); i++) {
if (model::cells[i] == c) {
model::cell_map[id] = i;
break;
}
}
@ -830,17 +836,17 @@ extern "C" {
void extend_cells_c(int32_t n)
{
cells.reserve(cells.size() + n);
model::cells.reserve(model::cells.size() + n);
for (int32_t i = 0; i < n; i++) {
cells.push_back(new CSGCell());
model::cells.push_back(new CSGCell());
}
n_cells = cells.size();
model::n_cells = model::cells.size();
}
int32_t universe_id(int i_univ) {return universes[i_univ]->id_;}
int32_t universe_id(int i_univ) {return model::universes[i_univ]->id_;}
void universes_to_hdf5(hid_t universes_group)
{for (Universe* u : universes) u->to_hdf5(universes_group);}
{for (Universe* u : model::universes) u->to_hdf5(universes_group);}
}

2
src/cmfd_execute.cpp
View file

@ -24,7 +24,7 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies,
openmc_source_bank(&source_bank, &n);
// Get source counts in each mesh bin / energy bin
auto& m = meshes.at(settings::index_cmfd_mesh);
auto& m = model::meshes.at(settings::index_cmfd_mesh);
xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
// Copy data from the xarray into the source counts array

5
src/constants.F90
View file

@ -15,7 +15,7 @@ module constants
VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE]
! HDF5 data format
integer, parameter :: HDF5_VERSION(2) = [1, 0]
integer, parameter :: HDF5_VERSION(2) = [2, 0]
! WMP data format
integer, parameter :: WMP_VERSION(2) = [1, 1]
@ -238,7 +238,8 @@ module constants
LIBRARY_NEUTRON = 1, &
LIBRARY_THERMAL = 2, &
LIBRARY_PHOTON = 3, &
LIBRARY_MULTIGROUP = 4
LIBRARY_MULTIGROUP = 4, &
LIBRARY_WMP = 5
! Probability table parameters
integer, parameter :: &

226
src/cross_sections.cpp Normal file
View file

@ -0,0 +1,226 @@
#include "openmc/cross_sections.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
#include "pugixml.hpp"
#include <cstdlib> // for getenv
namespace openmc {
//==============================================================================
// Global variable declarations
//==============================================================================
namespace data {
std::vector<Library> libraries;
std::map<LibraryKey, std::size_t> library_map;
}
//==============================================================================
// Library methods
//==============================================================================
Library::Library(pugi::xml_node node, const std::string& directory)
{
// Get type of library
if (check_for_node(node, "type")) {
auto type = get_node_value(node, "type");
if (type == "neutron") {
type_ = Type::neutron;
} else if (type == "thermal") {
type_ = Type::thermal;
} else if (type == "photon") {
type_ = Type::photon;
} else if (type == "wmp") {
type_ = Type::wmp;
} else {
fatal_error("Unrecognized library type: " + type);
}
} else {
fatal_error("Missing library type");
}
// Get list of materials
if (check_for_node(node, "materials")) {
materials_ = get_node_array<std::string>(node, "materials");
}
// determine path of cross section table
if (!check_for_node(node, "path")) {
fatal_error("Missing library path");
}
std::string path = get_node_value(node, "path");
if (starts_with(path, "/")) {
path_ = path;
} else if (ends_with(directory, "/")) {
path_ = directory + path;
} else {
path_ = directory + "/" + path;
}
if (!file_exists(path_)) {
warning("Cross section library " + path_ + " does not exist.");
}
}
//==============================================================================
// Non-member functions
//==============================================================================
extern "C" void read_mg_cross_sections_header();
void read_cross_sections_xml()
{
// Check if materials.xml exists
std::string filename = settings::path_input + "materials.xml";
if (!file_exists(filename)) {
fatal_error("Material XML file '" + filename + "' does not exist.");
}
// Parse materials.xml file
pugi::xml_document doc;
doc.load_file(filename.c_str());
auto root = doc.document_element();
// Find cross_sections.xml file -- the first place to look is the
// materials.xml file. If no file is found there, then we check the
// OPENMC_CROSS_SECTIONS environment variable
if (!check_for_node(root, "cross_sections")) {
// No cross_sections.xml file specified in settings.xml, check
// environment variable
if (settings::run_CE) {
char* envvar = std::getenv("OPENMC_CROSS_SECTIONS");
if (!envvar) {
fatal_error("No cross_sections.xml file was specified in "
"materials.xml or in the OPENMC_CROSS_SECTIONS"
" environment variable. OpenMC needs such a file to identify "
"where to find data libraries. Please consult the"
" user's guide at https://openmc.readthedocs.io for "
"information on how to set up data libraries.");
}
settings::path_cross_sections = envvar;
} else {
char* envvar = std::getenv("OPENMC_MG_CROSS_SECTIONS");
if (!envvar) {
fatal_error("No mgxs.h5 file was specified in "
"materials.xml or in the OPENMC_MG_CROSS_SECTIONS environment "
"variable. OpenMC needs such a file to identify where to "
"find MG cross section libraries. Please consult the user's "
"guide at http://openmc.readthedocs.io for information on "
"how to set up MG cross section libraries.");
}
settings::path_cross_sections = envvar;
}
} else {
settings::path_cross_sections = get_node_value(root, "cross_sections");
}
// Now that the cross_sections.xml or mgxs.h5 has been located, read it in
if (settings::run_CE) {
read_ce_cross_sections_xml();
} else {
read_mg_cross_sections_header();
}
// Establish mapping between (type, material) and index in libraries
int i = 0;
for (const auto& lib : data::libraries) {
for (const auto& name : lib.materials_) {
std::string lower_name = name;
to_lower(lower_name);
LibraryKey key {lib.type_, lower_name};
data::library_map.insert({key, i});
}
++i;
}
// Check that 0K nuclides are listed in the cross_sections.xml file
for (const auto& name : settings::res_scat_nuclides) {
std::string lower_name = name;
to_lower(lower_name);
LibraryKey key {Library::Type::neutron, lower_name};
if (data::library_map.find(key) == data::library_map.end()) {
fatal_error("Could not find resonant scatterer " +
name + " in cross_sections.xml file!");
}
}
}
void read_ce_cross_sections_xml()
{
// Check if cross_sections.xml exists
const auto& filename = settings::path_cross_sections;
if (!file_exists(filename)) {
// Could not find cross_sections.xml file
fatal_error("Cross sections XML file '" + filename +
"' does not exist.");
}
write_message("Reading cross sections XML file...", 5);
// Parse cross_sections.xml file
pugi::xml_document doc;
auto result = doc.load_file(filename.c_str());
if (!result) {
fatal_error("Error processing cross_sections.xml file.");
}
auto root = doc.document_element();
std::string directory;
if (check_for_node(root, "directory")) {
// Copy directory information if present
directory = get_node_value(root, "directory");
} else {
// If no directory is listed in cross_sections.xml, by default select the
// directory in which the cross_sections.xml file resides
auto pos = filename.rfind("/");
directory = filename.substr(0, pos);
}
for (const auto& node_library : root.children("library")) {
data::libraries.emplace_back(node_library, directory);
}
// Make sure file was not empty
if (data::libraries.empty()) {
fatal_error("No cross section libraries present in cross_sections.xml file.");
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void library_clear() {
data::libraries.clear();
data::library_map.clear();
}
extern "C" const char* library_path(int type, const char* name) {
auto lib_type = static_cast<Library::Type>(type);
LibraryKey key {lib_type, name};
if (data::library_map.find(key) == data::library_map.end()) {
return nullptr;
} else {
auto idx = data::library_map[key];
return data::libraries[idx].path_.c_str();
}
}
extern "C" bool library_present(int type, const char* name) {
auto lib_type = static_cast<Library::Type>(type);
LibraryKey key {lib_type, name};
return data::library_map.find(key) != data::library_map.end();
}
} // namespace openmc

81
src/dagmc.cpp
View file

@ -1,7 +1,8 @@
#include "openmc/dagmc.h"
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/string_functions.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
#include "openmc/geometry.h"
@ -13,20 +14,24 @@
namespace openmc {
namespace model {
moab::DagMC* DAG;
} // namespace model
void load_dagmc_geometry()
{
if (!DAG) {
DAG = new moab::DagMC();
if (!model::DAG) {
model::DAG = new moab::DagMC();
}
int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC
moab::ErrorCode rval = DAG->load_file("dagmc.h5m");
moab::ErrorCode rval = model::DAG->load_file("dagmc.h5m");
MB_CHK_ERR_CONT(rval);
rval = DAG->init_OBBTree();
rval = model::DAG->init_OBBTree();
MB_CHK_ERR_CONT(rval);
std::vector<std::string> prop_keywords;
@ -34,48 +39,45 @@ void load_dagmc_geometry()
prop_keywords.push_back("boundary");
std::map<std::string, std::string> ph;
DAG->parse_properties(prop_keywords, ph, ":");
model::DAG->parse_properties(prop_keywords, ph, ":");
MB_CHK_ERR_CONT(rval);
// initialize cell objects
n_cells = DAG->num_entities(3);
model::n_cells = model::DAG->num_entities(3);
// Allocate the cell overlap count if necessary.
if (settings::check_overlaps) overlap_check_count.resize(n_cells, 0);
for (int i = 0; i < n_cells; i++) {
moab::EntityHandle vol_handle = DAG->entity_by_index(3, i+1);
for (int i = 0; i < model::n_cells; i++) {
moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1);
// set cell ids using global IDs
DAGCell* c = new DAGCell();
c->id_ = DAG->id_by_index(3, i+1);
c->dagmc_ptr_ = DAG;
c->id_ = model::DAG->id_by_index(3, i+1);
c->dagmc_ptr_ = model::DAG;
c->universe_ = dagmc_univ_id; // set to zero for now
c->fill_ = C_NONE; // no fill, single universe
cells.push_back(c);
cell_map[c->id_] = i;
model::cells.push_back(c);
model::cell_map[c->id_] = i;
// Populate the Universe vector and dict
auto it = universe_map.find(dagmc_univ_id);
if (it == universe_map.end()) {
universes.push_back(new Universe());
universes.back()-> id_ = dagmc_univ_id;
universes.back()->cells_.push_back(i);
universe_map[dagmc_univ_id] = universes.size() - 1;
auto it = model::universe_map.find(dagmc_univ_id);
if (it == model::universe_map.end()) {
model::universes.push_back(new Universe());
model::universes.back()-> id_ = dagmc_univ_id;
model::universes.back()->cells_.push_back(i);
model::universe_map[dagmc_univ_id] = model::universes.size() - 1;
} else {
universes[it->second]->cells_.push_back(i);
model::universes[it->second]->cells_.push_back(i);
}
if (DAG->is_implicit_complement(vol_handle)) {
if (model::DAG->is_implicit_complement(vol_handle)) {
// assuming implicit complement is void for now
c->material_.push_back(MATERIAL_VOID);
continue;
}
if (DAG->has_prop(vol_handle, "mat")){
if (model::DAG->has_prop(vol_handle, "mat")){
std::string mat_value;
rval = DAG->prop_value(vol_handle, "mat", mat_value);
rval = model::DAG->prop_value(vol_handle, "mat", mat_value);
MB_CHK_ERR_CONT(rval);
to_lower(mat_value);
@ -91,21 +93,26 @@ void load_dagmc_geometry()
}
}
// Allocate the cell overlap count if necessary.
if (settings::check_overlaps) {
model::overlap_check_count.resize(model::cells.size(), 0);
}
// initialize surface objects
n_surfaces = DAG->num_entities(2);
surfaces.resize(n_surfaces);
int n_surfaces = model::DAG->num_entities(2);
model::surfaces.resize(n_surfaces);
for (int i = 0; i < n_surfaces; i++) {
moab::EntityHandle surf_handle = DAG->entity_by_index(2, i+1);
moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1);
// set cell ids using global IDs
DAGSurface* s = new DAGSurface();
s->id_ = DAG->id_by_index(2, i+1);
s->dagmc_ptr_ = DAG;
s->id_ = model::DAG->id_by_index(2, i+1);
s->dagmc_ptr_ = model::DAG;
if (DAG->has_prop(surf_handle, "boundary")) {
if (model::DAG->has_prop(surf_handle, "boundary")) {
std::string bc_value;
rval = DAG->prop_value(surf_handle, "boundary", bc_value);
rval = model::DAG->prop_value(surf_handle, "boundary", bc_value);
MB_CHK_ERR_CONT(rval);
to_lower(bc_value);
@ -128,8 +135,8 @@ void load_dagmc_geometry()
}
// add to global array and map
surfaces[i] = s;
surface_map[s->id_] = s->id_;
model::surfaces[i] = s;
model::surface_map[s->id_] = s->id_;
}
return;
@ -137,7 +144,7 @@ void load_dagmc_geometry()
void free_memory_dagmc()
{
delete DAG;
delete model::DAG;
}
}

78
src/eigenvalue.cpp
View file

@ -30,11 +30,15 @@ namespace openmc {
// Global variables
//==============================================================================
namespace simulation {
double keff_generation;
std::array<double, 2> k_sum;
std::vector<double> entropy;
xt::xtensor<double, 1> source_frac;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
@ -44,15 +48,15 @@ void calculate_generation_keff()
auto gt = global_tallies();
// Get keff for this generation by subtracting off the starting value
keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - keff_generation;
simulation::keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - simulation::keff_generation;
double keff_reduced;
#ifdef OPENMC_MPI
// Combine values across all processors
MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE,
MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE,
MPI_SUM, mpi::intracomm);
#else
keff_reduced = keff_generation;
keff_reduced = simulation::keff_generation;
#endif
// Normalize single batch estimate of k
@ -63,7 +67,7 @@ void calculate_generation_keff()
void synchronize_bank()
{
time_bank.start();
simulation::time_bank.start();
// Get pointers to source/fission bank
Bank* source_bank;
@ -83,21 +87,21 @@ void synchronize_bank()
#ifdef OPENMC_MPI
int64_t start = 0;
MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
// While we would expect the value of start on rank 0 to be 0, the MPI
// standard says that the receive buffer on rank 0 is undefined and not
// significant
if (mpi::rank == 0) start = 0;
int64_t finish = start + n_bank;
int64_t finish = start + simulation::n_bank;
int64_t total = finish;
MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
#else
int64_t start = 0;
int64_t finish = n_bank;
int64_t total = n_bank;
int64_t finish = simulation::n_bank;
int64_t total = simulation::n_bank;
#endif
// If there are not that many particles per generation, it's possible that no
@ -105,7 +109,7 @@ void synchronize_bank()
// extra logic to treat this circumstance, we really want to ensure the user
// runs enough particles to avoid this in the first place.
if (n_bank == 0) {
if (simulation::n_bank == 0) {
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
@ -127,7 +131,7 @@ void synchronize_bank()
}
double p_sample = static_cast<double>(sites_needed) / total;
time_bank_sample.start();
simulation::time_bank_sample.start();
// ==========================================================================
// SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
@ -136,7 +140,7 @@ void synchronize_bank()
int64_t index_temp = 0;
std::vector<Bank> temp_sites(3*simulation::work);
for (int64_t i = 0; i < n_bank; ++i) {
for (int64_t i = 0; i < simulation::n_bank; ++i) {
// If there are less than n_particles particles banked, automatically add
// int(n_particles/total) sites to temp_sites. For example, if you need
// 1000 and 300 were banked, this would add 3 source sites per banked site
@ -191,7 +195,7 @@ void synchronize_bank()
// fission bank
sites_needed = settings::n_particles - finish;
for (int i = 0; i < sites_needed; ++i) {
temp_sites[index_temp] = fission_bank[n_bank - sites_needed + i];
temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i];
++index_temp;
}
}
@ -200,8 +204,8 @@ void synchronize_bank()
finish = simulation::work_index[mpi::rank + 1];
}
time_bank_sample.stop();
time_bank_sendrecv.start();
simulation::time_bank_sample.stop();
simulation::time_bank_sendrecv.start();
#ifdef OPENMC_MPI
// ==========================================================================
@ -299,8 +303,8 @@ void synchronize_bank()
std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank);
#endif
time_bank_sendrecv.stop();
time_bank.stop();
simulation::time_bank_sendrecv.stop();
simulation::time_bank.stop();
}
void calculate_average_keff()
@ -320,11 +324,11 @@ void calculate_average_keff()
simulation::keff = simulation::k_generation[i];
} else {
// Sample mean of keff
k_sum[0] += simulation::k_generation[i];
k_sum[1] += std::pow(simulation::k_generation[i], 2);
simulation::k_sum[0] += simulation::k_generation[i];
simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
// Determine mean
simulation::keff = k_sum[0] / n;
simulation::keff = simulation::k_sum[0] / n;
if (n > 1) {
double t_value;
@ -337,7 +341,7 @@ void calculate_average_keff()
}
// Standard deviation of the sample mean of k
simulation::keff_std = t_value * std::sqrt((k_sum[1]/n -
simulation::keff_std = t_value * std::sqrt((simulation::k_sum[1]/n -
std::pow(simulation::keff, 2)) / (n - 1));
}
}
@ -493,7 +497,7 @@ int openmc_get_keff(double* k_combined)
void shannon_entropy()
{
// Get pointer to entropy mesh
auto& m = meshes[settings::index_entropy_mesh];
auto& m = model::meshes[settings::index_entropy_mesh];
// Get pointer to fission bank
Bank* fission_bank;
@ -503,7 +507,7 @@ void shannon_entropy()
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(
n_bank, fission_bank, 0, nullptr, &sites_outside);
simulation::n_bank, fission_bank, 0, nullptr, &sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
@ -523,20 +527,20 @@ void shannon_entropy()
}
// Add value to vector
entropy.push_back(H);
simulation::entropy.push_back(H);
}
}
void ufs_count_sites()
{
auto &m = meshes[settings::index_ufs_mesh];
auto &m = model::meshes[settings::index_ufs_mesh];
if (simulation::current_batch == 1 && simulation::current_gen == 1) {
// On the first generation, just assume that the source is already evenly
// distributed so that effectively the production of fission sites is not
// biased
auto s = xt::view(source_frac, xt::all());
auto s = xt::view(simulation::source_frac, xt::all());
s = m->volume_frac_;
} else {
@ -547,7 +551,7 @@ void ufs_count_sites()
// count number of source sites in each ufs mesh cell
bool sites_outside;
source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
&sites_outside);
// Check for sites outside of the mesh
@ -558,12 +562,12 @@ void ufs_count_sites()
#ifdef OPENMC_MPI
// Send source fraction to all processors
int n_bins = xt::prod(m->shape_)();
MPI_Bcast(source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
MPI_Bcast(simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
#endif
// Normalize to total weight to get fraction of source in each cell
double total = xt::sum(source_frac)();
source_frac /= total;
double total = xt::sum(simulation::source_frac)();
simulation::source_frac /= total;
// Since the total starting weight is not equal to n_particles, we need to
// renormalize the weight of the source sites
@ -575,7 +579,7 @@ void ufs_count_sites()
double ufs_get_weight(const Particle* p)
{
auto& m = meshes[settings::index_ufs_mesh];
auto& m = model::meshes[settings::index_ufs_mesh];
// Determine indices on ufs mesh for current location
// TODO: off by one
@ -585,8 +589,8 @@ double ufs_get_weight(const Particle* p)
fatal_error("Source site outside UFS mesh!");
}
if (source_frac(mesh_bin) != 0.0) {
return m->volume_frac_ / source_frac(mesh_bin);
if (simulation::source_frac(mesh_bin) != 0.0) {
return m->volume_frac_ / simulation::source_frac(mesh_bin);
} else {
return 1.0;
}
@ -598,7 +602,7 @@ extern "C" void write_eigenvalue_hdf5(hid_t group)
write_dataset(group, "generations_per_batch", settings::gen_per_batch);
write_dataset(group, "k_generation", simulation::k_generation);
if (settings::entropy_on) {
write_dataset(group, "entropy", entropy);
write_dataset(group, "entropy", simulation::entropy);
}
write_dataset(group, "k_col_abs", simulation::k_col_abs);
write_dataset(group, "k_col_tra", simulation::k_col_tra);
@ -615,7 +619,7 @@ extern "C" void read_eigenvalue_hdf5(hid_t group)
simulation::k_generation.resize(n);
read_dataset(group, "k_generation", simulation::k_generation);
if (settings::entropy_on) {
read_dataset(group, "entropy", entropy);
read_dataset(group, "entropy", simulation::entropy);
}
read_dataset(group, "k_col_abs", simulation::k_col_abs);
read_dataset(group, "k_col_tra", simulation::k_col_tra);
@ -628,14 +632,14 @@ extern "C" void read_eigenvalue_hdf5(hid_t group)
extern "C" double entropy_c(int i)
{
return entropy.at(i - 1);
return simulation::entropy.at(i - 1);
}
extern "C" void entropy_clear()
{
entropy.clear();
simulation::entropy.clear();
}
extern "C" void k_sum_reset() { k_sum.fill(0.0); }
extern "C" void k_sum_reset() { simulation::k_sum.fill(0.0); }
} // namespace openmc

22
src/endf.F90
View file

@ -229,7 +229,7 @@ contains
integer, intent(in) :: MT
logical :: dis
if (MT >= N_GAMMA .and. MT <= N_DA) then
if (MT >= N_DISAPPEAR .and. MT <= N_DA) then
dis = .true.
elseif (MT >= N_P0 .and. MT <= N_AC) then
dis = .true.
@ -242,25 +242,27 @@ contains
end function is_disappearance
!===============================================================================
! IS_SCATTER determines if a given MT number is that of a scattering event
! IS_INELASTIC_SCATTER determines if a given MT number is that of an inelastic
! scattering event
!===============================================================================
function is_scatter(MT) result(scatter_event)
function is_inelastic_scatter(MT) result(retval)
integer, intent(in) :: MT
logical :: scatter_event
logical :: retval
if (MT < 100) then
if (MT == N_FISSION .or. MT == N_F .or. MT == N_NF .or. MT == N_2NF &
.or. MT == N_3NF) then
scatter_event = .false.
if (is_fission(MT)) then
retval = .false.
else
scatter_event = .true.
retval = (MT >= MISC .and. MT /= 27)
end if
elseif (MT <= 200) then
retval = (.not. is_disappearance(MT))
else
scatter_event = .false.
retval = .false.
end if
end function is_scatter
end function
end module endf

8
src/finalize.cpp
View file

@ -75,10 +75,10 @@ int openmc_finalize()
simulation::satisfy_triggers = false;
simulation::total_gen = 0;
energy_max = {INFTY, INFTY};
energy_min = {0.0, 0.0};
data::energy_max = {INFTY, INFTY};
data::energy_min = {0.0, 0.0};
n_tallies = 0;
openmc_root_universe = -1;
model::root_universe = -1;
openmc_set_seed(DEFAULT_SEED);
// Deallocate arrays
@ -113,7 +113,7 @@ int openmc_reset()
simulation::k_col_abs = 0.0;
simulation::k_col_tra = 0.0;
simulation::k_abs_tra = 0.0;
k_sum = {0.0, 0.0};
simulation::k_sum = {0.0, 0.0};
// Reset timers
reset_timers();

53
src/geometry.cpp
View file

@ -14,8 +14,21 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace model {
int root_universe {-1};
std::vector<int64_t> overlap_check_count;
} // namespace model
//==============================================================================
// Non-member functions
//==============================================================================
extern "C" bool
@ -24,21 +37,21 @@ check_cell_overlap(Particle* p) {
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *universes[p->coord[j].universe];
Universe& univ = *model::universes[p->coord[j].universe];
int n = univ.cells_.size();
// Loop through each cell on this level
for (auto index_cell : univ.cells_) {
Cell& c = *cells[index_cell];
Cell& c = *model::cells[index_cell];
if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) {
if (index_cell != p->coord[j].cell) {
std::stringstream err_msg;
err_msg << "Overlapping cells detected: " << c.id_ << ", "
<< cells[p->coord[j].cell]->id_ << " on universe "
<< model::cells[p->coord[j].cell]->id_ << " on universe "
<< univ.id_;
fatal_error(err_msg);
}
++overlap_check_count[index_cell];
++model::overlap_check_count[index_cell];
}
}
}
@ -57,8 +70,8 @@ find_cell(Particle* p, int search_surf) {
// Determine universe (if not yet set, use root universe)
int i_universe = p->coord[p->n_coord-1].universe;
if (i_universe == C_NONE) {
p->coord[p->n_coord-1].universe = openmc_root_universe;
i_universe = openmc_root_universe;
p->coord[p->n_coord-1].universe = model::root_universe;
i_universe = model::root_universe;
}
// If a surface was indicated, only search cells from the neighbor list of
@ -66,12 +79,12 @@ find_cell(Particle* p, int search_surf) {
// the positive or negative side of the surface should be searched.
const std::vector<int>* search_cells;
if (search_surf > 0) {
search_cells = &surfaces[search_surf-1]->neighbor_pos_;
search_cells = &model::surfaces[search_surf-1]->neighbor_pos_;
} else if (search_surf < 0) {
search_cells = &surfaces[-search_surf-1]->neighbor_neg_;
search_cells = &model::surfaces[-search_surf-1]->neighbor_neg_;
} else {
// No surface was indicated, search all cells in the universe.
search_cells = &universes[i_universe]->cells_;
search_cells = &model::universes[i_universe]->cells_;
}
// Find which cell of this universe the particle is in.
@ -81,17 +94,17 @@ find_cell(Particle* p, int search_surf) {
i_cell = (*search_cells)[i];
// Make sure the search cell is in the same universe.
if (cells[i_cell]->universe_ != i_universe) continue;
if (model::cells[i_cell]->universe_ != i_universe) continue;
Position r {p->coord[p->n_coord-1].xyz};
Direction u {p->coord[p->n_coord-1].uvw};
int32_t surf = p->surface;
if (cells[i_cell]->contains(r, u, surf)) {
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord[p->n_coord-1].cell = i_cell;
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
msg << " Entering cell " << cells[i_cell]->id_;
msg << " Entering cell " << model::cells[i_cell]->id_;
write_message(msg, 1);
}
found = true;
@ -100,7 +113,7 @@ find_cell(Particle* p, int search_surf) {
}
if (found) {
Cell& c {*cells[i_cell]};
Cell& c {*model::cells[i_cell]};
if (c.type_ == FILL_MATERIAL) {
//=======================================================================
//! Found a material cell which means this is the lowest coord level.
@ -109,11 +122,11 @@ find_cell(Particle* p, int search_surf) {
if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) {
int offset = 0;
for (int i = 0; i < p->n_coord; i++) {
Cell& c_i {*cells[p->coord[i].cell]};
Cell& c_i {*model::cells[p->coord[i].cell]};
if (c_i.type_ == FILL_UNIVERSE) {
offset += c_i.offset_[c.distribcell_index_];
} else if (c_i.type_ == FILL_LATTICE) {
Lattice& lat {*lattices[p->coord[i+1].lattice-1]};
Lattice& lat {*model::lattices[p->coord[i+1].lattice-1]};
int i_xyz[3] {p->coord[i+1].lattice_x,
p->coord[i+1].lattice_y,
p->coord[i+1].lattice_z};
@ -198,7 +211,7 @@ find_cell(Particle* p, int search_surf) {
//========================================================================
//! Found a lower lattice, update this coord level then search the next.
Lattice& lat {*lattices[c.fill_]};
Lattice& lat {*model::lattices[c.fill_]};
// Determine lattice indices.
Position r {p->coord[p->n_coord-1].xyz};
@ -251,7 +264,7 @@ find_cell(Particle* p, int search_surf) {
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3])
{
Lattice& lat {*lattices[p->coord[p->n_coord-1].lattice-1]};
Lattice& lat {*model::lattices[p->coord[p->n_coord-1].lattice-1]};
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
@ -326,7 +339,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
for (int i = 0; i < p->n_coord; i++) {
Position r {p->coord[i].xyz};
Direction u {p->coord[i].uvw};
Cell& c {*cells[p->coord[i].cell]};
Cell& c {*model::cells[p->coord[i].cell]};
// Find the oncoming surface in this cell and the distance to it.
auto surface_distance = c.distance(r, u, p->surface);
@ -335,7 +348,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
// Find the distance to the next lattice tile crossing.
if (p->coord[i].lattice != F90_NONE) {
Lattice& lat {*lattices[p->coord[i].lattice-1]};
Lattice& lat {*model::lattices[p->coord[i].lattice-1]};
std::array<int, 3> i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y,
p->coord[i].lattice_z};
//TODO: refactor so both lattice use the same position argument (which
@ -377,7 +390,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
*surface_crossed = level_surf_cross;
} else {
Position r_hit = r + d_surf * u;
Surface& surf {*surfaces[std::abs(level_surf_cross)-1]};
Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]};
Direction norm = surf.normal(r_hit);
if (u.dot(norm) > 0) {
*surface_crossed = std::abs(level_surf_cross);

124
src/geometry_aux.cpp
View file

@ -24,15 +24,15 @@ void
adjust_indices()
{
// Adjust material/fill idices.
for (Cell* c : cells) {
for (Cell* c : model::cells) {
if (c->fill_ != C_NONE) {
int32_t id = c->fill_;
auto search_univ = universe_map.find(id);
auto search_lat = lattice_map.find(id);
if (search_univ != universe_map.end()) {
auto search_univ = model::universe_map.find(id);
auto search_lat = model::lattice_map.find(id);
if (search_univ != model::universe_map.end()) {
c->type_ = FILL_UNIVERSE;
c->fill_ = search_univ->second;
} else if (search_lat != lattice_map.end()) {
} else if (search_lat != model::lattice_map.end()) {
c->type_ = FILL_LATTICE;
c->fill_ = search_lat->second;
} else {
@ -46,8 +46,8 @@ adjust_indices()
for (auto it = c->material_.begin(); it != c->material_.end(); it++) {
int32_t mid = *it;
if (mid != MATERIAL_VOID) {
auto search = material_map.find(mid);
if (search != material_map.end()) {
auto search = model::material_map.find(mid);
if (search != model::material_map.end()) {
*it = search->second;
} else {
std::stringstream err_msg;
@ -61,9 +61,9 @@ adjust_indices()
}
// Change cell.universe values from IDs to indices.
for (Cell* c : cells) {
auto search = universe_map.find(c->universe_);
if (search != universe_map.end()) {
for (Cell* c : model::cells) {
auto search = model::universe_map.find(c->universe_);
if (search != model::universe_map.end()) {
c->universe_ = search->second;
} else {
std::stringstream err_msg;
@ -74,7 +74,7 @@ adjust_indices()
}
// Change all lattice universe values from IDs to indices.
for (Lattice* l : lattices) {
for (Lattice* l : model::lattices) {
l->adjust_indices();
}
}
@ -84,7 +84,7 @@ adjust_indices()
void
assign_temperatures()
{
for (Cell* c : cells) {
for (Cell* c : model::cells) {
// Ignore non-material cells and cells with defined temperature.
if (c->material_.size() == 0) continue;
if (c->sqrtkT_.size() > 0) continue;
@ -96,9 +96,9 @@ assign_temperatures()
c->sqrtkT_.push_back(0);
} else {
if (materials[i_mat]->temperature_ >= 0) {
if (model::materials[i_mat]->temperature_ >= 0) {
// This material has a default temperature; use that value.
auto T = materials[i_mat]->temperature_;
auto T = model::materials[i_mat]->temperature_;
c->sqrtkT_.push_back(std::sqrt(K_BOLTZMANN * T));
} else {
// Use the global default temperature.
@ -117,12 +117,12 @@ find_root_universe()
{
// Find all the universes listed as a cell fill.
std::unordered_set<int32_t> fill_univ_ids;
for (Cell* c : cells) {
for (Cell* c : model::cells) {
fill_univ_ids.insert(c->fill_);
}
// Find all the universes contained in a lattice.
for (Lattice* lat : lattices) {
for (Lattice* lat : model::lattices) {
for (auto it = lat->begin(); it != lat->end(); ++it) {
fill_univ_ids.insert(*it);
}
@ -134,8 +134,8 @@ find_root_universe()
// Figure out which universe is not in the set. This is the root universe.
bool root_found {false};
int32_t root_univ;
for (int32_t i = 0; i < universes.size(); i++) {
auto search = fill_univ_ids.find(universes[i]->id_);
for (int32_t i = 0; i < model::universes.size(); i++) {
auto search = fill_univ_ids.find(model::universes[i]->id_);
if (search == fill_univ_ids.end()) {
if (root_found) {
fatal_error("Two or more universes are not used as fill universes, so "
@ -160,21 +160,21 @@ neighbor_lists()
{
write_message("Building neighboring cells lists for each surface...", 6);
for (int i = 0; i < cells.size(); i++) {
for (auto token : cells[i]->region_) {
for (int i = 0; i < model::cells.size(); i++) {
for (auto token : model::cells[i]->region_) {
// Skip operator tokens.
if (std::abs(token) >= OP_UNION) continue;
// This token is a surface index. Add the cell to the surface's list.
if (token > 0) {
surfaces[std::abs(token)-1]->neighbor_pos_.push_back(i);
model::surfaces[std::abs(token)-1]->neighbor_pos_.push_back(i);
} else {
surfaces[std::abs(token)-1]->neighbor_neg_.push_back(i);
model::surfaces[std::abs(token)-1]->neighbor_neg_.push_back(i);
}
}
}
for (Surface* surf : surfaces) {
for (Surface* surf : model::surfaces) {
surf->neighbor_pos_.shrink_to_fit();
surf->neighbor_neg_.shrink_to_fit();
}
@ -187,7 +187,7 @@ prepare_distribcell()
{
// Find all cells listed in a DistribcellFilter.
std::unordered_set<int32_t> distribcells;
for (auto& filt : tally_filters) {
for (auto& filt : model::tally_filters) {
auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
if (distrib_filt) {
distribcells.insert(distrib_filt->cell_);
@ -196,8 +196,8 @@ prepare_distribcell()
// Find all cells with distributed materials or temperatures. Make sure that
// the number of materials/temperatures matches the number of cell instances.
for (int i = 0; i < cells.size(); i++) {
Cell& c {*cells[i]};
for (int i = 0; i < model::cells.size(); i++) {
Cell& c {*model::cells[i]};
if (c.material_.size() > 1) {
if (c.material_.size() != c.n_instances_) {
@ -228,10 +228,10 @@ prepare_distribcell()
// unique distribcell array index.
int distribcell_index = 0;
std::vector<int32_t> target_univ_ids;
for (Universe* u : universes) {
for (Universe* u : model::universes) {
for (auto cell_indx : u->cells_) {
if (distribcells.find(cell_indx) != distribcells.end()) {
cells[cell_indx]->distribcell_index_ = distribcell_index;
model::cells[cell_indx]->distribcell_index_ = distribcell_index;
target_univ_ids.push_back(u->id_);
++distribcell_index;
}
@ -240,22 +240,22 @@ prepare_distribcell()
// Allocate the cell and lattice offset tables.
int n_maps = target_univ_ids.size();
for (Cell* c : cells) {
for (Cell* c : model::cells) {
if (c->type_ != FILL_MATERIAL) {
c->offset_.resize(n_maps, C_NONE);
}
}
for (Lattice* lat : lattices) {
for (Lattice* lat : model::lattices) {
lat->allocate_offset_table(n_maps);
}
// Fill the cell and lattice offset tables.
for (int map = 0; map < target_univ_ids.size(); map++) {
auto target_univ_id = target_univ_ids[map];
for (Universe* univ : universes) {
for (Universe* univ : model::universes) {
int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation.
for (int32_t cell_indx : univ->cells_) {
Cell& c = *cells[cell_indx];
Cell& c = *model::cells[cell_indx];
if (c.type_ == FILL_UNIVERSE) {
c.offset_[map] = offset;
@ -263,7 +263,7 @@ prepare_distribcell()
offset += count_universe_instances(search_univ, target_univ_id);
} else if (c.type_ == FILL_LATTICE) {
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
offset = lat.fill_offset_table(offset, target_univ_id, map);
}
}
@ -276,8 +276,8 @@ prepare_distribcell()
void
count_cell_instances(int32_t univ_indx)
{
for (int32_t cell_indx : universes[univ_indx]->cells_) {
Cell& c = *cells[cell_indx];
for (int32_t cell_indx : model::universes[univ_indx]->cells_) {
Cell& c = *model::cells[cell_indx];
++c.n_instances_;
if (c.type_ == FILL_UNIVERSE) {
@ -286,7 +286,7 @@ count_cell_instances(int32_t univ_indx)
} else if (c.type_ == FILL_LATTICE) {
// This cell contains a lattice. Recurse into the lattice universes.
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
for (auto it = lat.begin(); it != lat.end(); ++it) {
count_cell_instances(*it);
}
@ -300,20 +300,20 @@ int
count_universe_instances(int32_t search_univ, int32_t target_univ_id)
{
// If this is the target, it can't contain itself.
if (universes[search_univ]->id_ == target_univ_id) {
if (model::universes[search_univ]->id_ == target_univ_id) {
return 1;
}
int count {0};
for (int32_t cell_indx : universes[search_univ]->cells_) {
Cell& c = *cells[cell_indx];
for (int32_t cell_indx : model::universes[search_univ]->cells_) {
Cell& c = *model::cells[cell_indx];
if (c.type_ == FILL_UNIVERSE) {
int32_t next_univ = c.fill_;
count += count_universe_instances(next_univ, target_univ_id);
} else if (c.type_ == FILL_LATTICE) {
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
for (auto it = lat.begin(); it != lat.end(); ++it) {
int32_t next_univ = *it;
count += count_universe_instances(next_univ, target_univ_id);
@ -338,7 +338,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
// write to the path and return.
for (int32_t cell_indx : search_univ.cells_) {
if ((cell_indx == target_cell) && (offset == target_offset)) {
Cell& c = *cells[cell_indx];
Cell& c = *model::cells[cell_indx];
path << "c" << c.id_;
return path.str();
}
@ -350,7 +350,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
std::vector<std::int32_t>::const_reverse_iterator cell_it
{search_univ.cells_.crbegin()};
for (; cell_it != search_univ.cells_.crend(); ++cell_it) {
Cell& c = *cells[*cell_it];
Cell& c = *model::cells[*cell_it];
// Material cells don't contain other cells so ignore them.
if (c.type_ != FILL_MATERIAL) {
@ -358,7 +358,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
if (c.type_ == FILL_UNIVERSE) {
temp_offset = offset + c.offset_[map];
} else {
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
int32_t indx = lat.universes_.size()*map + lat.begin().indx_;
temp_offset = offset + lat.offsets_[indx];
}
@ -370,18 +370,18 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
}
// Add the cell to the path string.
Cell& c = *cells[*cell_it];
Cell& c = *model::cells[*cell_it];
path << "c" << c.id_ << "->";
if (c.type_ == FILL_UNIVERSE) {
// Recurse into the fill cell.
offset += c.offset_[map];
path << distribcell_path_inner(target_cell, map, target_offset,
*universes[c.fill_], offset);
*model::universes[c.fill_], offset);
return path.str();
} else {
// Recurse into the lattice cell.
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
path << "l" << lat.id_;
for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) {
int32_t indx = lat.universes_.size()*map + it.indx_;
@ -390,7 +390,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
offset = temp_offset;
path << "(" << lat.index_to_string(it.indx_) << ")->";
path << distribcell_path_inner(target_cell, map, target_offset,
*universes[*it], offset);
*model::universes[*it], offset);
return path.str();
}
}
@ -401,7 +401,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
std::string
distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset)
{
auto& root_univ = *universes[openmc_root_universe];
auto& root_univ = *model::universes[model::root_universe];
return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0);
}
@ -412,13 +412,13 @@ maximum_levels(int32_t univ)
{
int levels_below {0};
for (int32_t cell_indx : universes[univ]->cells_) {
Cell& c = *cells[cell_indx];
for (int32_t cell_indx : model::universes[univ]->cells_) {
Cell& c = *model::cells[cell_indx];
if (c.type_ == FILL_UNIVERSE) {
int32_t next_univ = c.fill_;
levels_below = std::max(levels_below, maximum_levels(next_univ));
} else if (c.type_ == FILL_LATTICE) {
Lattice& lat = *lattices[c.fill_];
Lattice& lat = *model::lattices[c.fill_];
for (auto it = lat.begin(); it != lat.end(); ++it) {
int32_t next_univ = *it;
levels_below = std::max(levels_below, maximum_levels(next_univ));
@ -435,20 +435,20 @@ maximum_levels(int32_t univ)
void
free_memory_geometry_c()
{
for (Cell* c : cells) {delete c;}
cells.clear();
cell_map.clear();
n_cells = 0;
for (Cell* c : model::cells) {delete c;}
model::cells.clear();
model::cell_map.clear();
model::n_cells = 0;
for (Universe* u : universes) {delete u;}
universes.clear();
universe_map.clear();
for (Universe* u : model::universes) {delete u;}
model::universes.clear();
model::universe_map.clear();
for (Lattice* lat : lattices) {delete lat;}
lattices.clear();
lattice_map.clear();
for (Lattice* lat : model::lattices) {delete lat;}
model::lattices.clear();
model::lattice_map.clear();
overlap_check_count.clear();
model::overlap_check_count.clear();
}
} // namespace openmc

2
src/geometry_header.F90
View file

@ -143,7 +143,7 @@ module geometry_header
end type Cell
! array index of the root universe
integer(C_INT), bind(C, name='openmc_root_universe') :: root_universe = -1
integer(C_INT), bind(C) :: root_universe
integer(C_INT32_T), bind(C) :: n_cells ! # of cells
integer(C_INT32_T), bind(C) :: n_universes ! # of universes

2
src/hdf5_interface.cpp
View file

@ -362,7 +362,7 @@ member_names(hid_t group_id, H5O_type_t type)
char buffer[size];
H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i,
buffer, size, H5P_DEFAULT);
names.emplace_back(&buffer[0], size);
names.emplace_back(&buffer[0]);
}
return names;
}

26
src/initialize.F90
View file

@ -8,23 +8,23 @@ module initialize
implicit none
interface
function openmc_path_input() result(ptr) bind(C)
function path_input_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_output() result(ptr) bind(C)
function path_output_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_particle_restart() result(ptr) bind(C)
function path_particle_restart_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_statepoint() result(ptr) bind(C)
function path_statepoint_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_sourcepoint() result(ptr) bind(C)
function path_sourcepoint_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
@ -49,22 +49,22 @@ contains
end function is_null
end interface
if (.not. is_null(openmc_path_input())) then
call c_f_pointer(openmc_path_input(), string, [255])
if (.not. is_null(path_input_c())) then
call c_f_pointer(path_input_c(), string, [255])
path_input = to_f_string(string)
else
path_input = ''
end if
if (.not. is_null(openmc_path_statepoint())) then
call c_f_pointer(openmc_path_statepoint(), string, [255])
if (.not. is_null(path_statepoint_c())) then
call c_f_pointer(path_statepoint_c(), string, [255])
path_state_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_sourcepoint())) then
call c_f_pointer(openmc_path_sourcepoint(), string, [255])
if (.not. is_null(path_sourcepoint_c())) then
call c_f_pointer(path_sourcepoint_c(), string, [255])
path_source_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_particle_restart())) then
call c_f_pointer(openmc_path_particle_restart(), string, [255])
if (.not. is_null(path_particle_restart_c())) then
call c_f_pointer(path_particle_restart_c(), string, [255])
path_particle_restart = to_f_string(string)
end if
end subroutine read_command_line

8
src/initialize.cpp
View file

@ -54,8 +54,8 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
if (err) return err;
// Start total and initialization timer
time_total.start();
time_initialize.start();
simulation::time_total.start();
simulation::time_initialize.start();
#ifdef _OPENMP
// If OMP_SCHEDULE is not set, default to a static schedule
@ -79,7 +79,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
if (settings::particle_restart_run) settings::run_mode = RUN_MODE_PARTICLE;
// Stop initialization timer
time_initialize.stop();
simulation::time_initialize.stop();
return 0;
}
@ -210,7 +210,7 @@ parse_command_line(int argc, char* argv[])
}
omp_set_num_threads(simulation::n_threads);
#else
if (openmc_master)
if (mpi::master)
warning("Ignoring number of threads specified on command line.");
#endif

370
src/input_xml.F90
View file

@ -72,6 +72,9 @@ module input_xml
type(C_PTR) :: node_ptr
end subroutine read_cells
subroutine read_cross_sections_xml() bind(C)
end subroutine
subroutine read_lattices(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
@ -123,15 +126,14 @@ contains
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
real(8), allocatable :: material_temps(:)
call read_settings_xml()
call read_cross_sections_xml()
call read_materials_xml(material_temps)
call read_materials_xml()
call read_geometry_xml()
! Set up neighbor lists, convert user IDs -> indices, assign temperatures
call finalize_geometry(material_temps, nuc_temps, sab_temps)
call finalize_geometry(nuc_temps, sab_temps)
if (run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
@ -172,8 +174,7 @@ contains
end subroutine read_input_xml
subroutine finalize_geometry(material_temps, nuc_temps, sab_temps)
real(8), intent(in) :: material_temps(:)
subroutine finalize_geometry(nuc_temps, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
@ -215,27 +216,12 @@ contains
integer :: i
integer :: n
type(XMLNode) :: root
type(XMLNode) :: node_res_scat
type(XMLNode) :: node_vol
type(XMLNode), allocatable :: node_vol_list(:)
! Get proper XMLNode type given pointer
root % ptr = root_ptr
! Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) then
node_res_scat = root % child("resonance_scattering")
! Get nuclides that resonance scattering should be applied to
if (check_for_node(node_res_scat, "nuclides")) then
n = node_word_count(node_res_scat, "nuclides")
allocate(res_scat_nuclides(n))
if (n > 0) then
call get_node_array(node_res_scat, "nuclides", res_scat_nuclides)
end if
end if
end if
call get_node_list(root, "volume_calc", node_vol_list)
n = size(node_vol_list)
allocate(volume_calcs(n))
@ -361,8 +347,8 @@ contains
call read_surfaces(root % ptr)
! Allocate surfaces array
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
allocate(surfaces(surfaces_size()))
do i = 1, size(surfaces)
surfaces(i) % ptr = surface_pointer(i - 1);
if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
@ -519,9 +505,9 @@ contains
integer :: i
! Allocate surfaces array
allocate(surfaces(n_surfaces))
allocate(surfaces(surfaces_size()))
do i = 1, n_surfaces
do i = 1, size(surfaces)
surfaces(i) % ptr = surface_pointer(i - 1);
end do
@ -547,128 +533,12 @@ contains
end do
end subroutine allocate_cells
!===============================================================================
! READ_MATERIAL_XML reads data from a materials.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_cross_sections_xml()
integer :: i, j
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable
character(MAX_LINE_LEN) :: filename
type(XMLDocument) :: doc
type(XMLNode) :: root
! Check if materials.xml exists
filename = trim(path_input) // "materials.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Material XML file '" // trim(filename) // "' does not &
&exist!")
end if
! Parse materials.xml file
call doc % load_file(filename)
root = doc % document_element()
! Find cross_sections.xml file -- the first place to look is the
! materials.xml file. If no file is found there, then we check the
! OPENMC_CROSS_SECTIONS environment variable
if (.not. check_for_node(root, "cross_sections")) then
! No cross_sections.xml file specified in settings.xml, check
! environment variable
if (run_CE) then
call get_environment_variable("OPENMC_CROSS_SECTIONS", env_variable)
if (len_trim(env_variable) == 0) then
call get_environment_variable("CROSS_SECTIONS", env_variable)
! FIXME: When deprecated option of setting the cross sections in
! settings.xml is removed, remove ".and. path_cross_sections == ''"
if (len_trim(env_variable) == 0 .and. path_cross_sections == '') then
call fatal_error("No cross_sections.xml file was specified in &
&materials.xml, settings.xml, or in the OPENMC_CROSS_SECTIONS&
& environment variable. OpenMC needs such a file to identify &
&where to find ACE cross section libraries. Please consult the&
& user's guide at http://openmc.readthedocs.io for &
&information on how to set up ACE cross section libraries.")
else
call warning("The CROSS_SECTIONS environment variable is &
&deprecated. Please update your environment to use &
&OPENMC_CROSS_SECTIONS instead.")
end if
end if
path_cross_sections = trim(env_variable)
else
call get_environment_variable("OPENMC_MG_CROSS_SECTIONS", env_variable)
! FIXME: When deprecated option of setting the mg cross sections in
! settings.xml is removed, remove ".and. path_cross_sections == ''"
if (len_trim(env_variable) == 0 .and. path_cross_sections == '') then
call fatal_error("No mgxs.h5 file was specified in &
&materials.xml or in the OPENMC_MG_CROSS_SECTIONS environment &
&variable. OpenMC needs such a file to identify where to &
&find MG cross section libraries. Please consult the user's &
&guide at http://openmc.readthedocs.io for information on &
&how to set up MG cross section libraries.")
else if (len_trim(env_variable) /= 0) then
path_cross_sections = trim(env_variable)
end if
end if
else
call get_node_value(root, "cross_sections", path_cross_sections)
end if
! Find the windowed multipole library
if (run_mode /= MODE_PLOTTING) then
if (.not. check_for_node(root, "multipole_library")) then
! No library location specified in materials.xml, check
! environment variable
call get_environment_variable("OPENMC_MULTIPOLE_LIBRARY", env_variable)
path_multipole = trim(env_variable)
else
call get_node_value(root, "multipole_library", path_multipole)
end if
if (.not. ends_with(path_multipole, "/")) &
path_multipole = trim(path_multipole) // "/"
end if
! Close materials XML file
call doc % clear()
! Now that the cross_sections.xml or mgxs.h5 has been located, read it in
if (run_CE) then
call read_ce_cross_sections_xml()
else
call read_mg_cross_sections_header()
end if
! Creating dictionary that maps the name of the material to the entry
do i = 1, size(libraries)
do j = 1, size(libraries(i) % materials)
call library_dict % set(to_lower(libraries(i) % materials(j)), i)
end do
end do
! Check that 0K nuclides are listed in the cross_sections.xml file
if (allocated(res_scat_nuclides)) then
do i = 1, size(res_scat_nuclides)
if (.not. library_dict % has(to_lower(res_scat_nuclides(i)))) then
call fatal_error("Could not find resonant scatterer " &
// trim(res_scat_nuclides(i)) // " in cross_sections.xml file!")
end if
end do
end if
end subroutine read_cross_sections_xml
subroutine read_materials_xml(material_temps)
real(8), allocatable, intent(out) :: material_temps(:)
subroutine read_materials_xml()
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: i_library ! index in libraries array
integer :: index_nuclide ! index in nuclides
integer :: index_element ! index in elements
integer :: index_sab ! index in sab_tables
@ -720,7 +590,6 @@ contains
! Allocate materials array
n_materials = size(node_mat_list)
allocate(materials(n_materials))
allocate(material_temps(n_materials))
! Initialize count for number of nuclides/S(a,b) tables
index_nuclide = 0
@ -745,13 +614,6 @@ contains
call get_node_value(node_mat, "name", mat % name)
end if
! Get material default temperature
if (check_for_node(node_mat, "temperature")) then
call get_node_value(node_mat, "temperature", material_temps(i))
else
material_temps(i) = -1.0
end if
! Get pointer to density element
if (check_for_node(node_mat, "density")) then
node_dens = node_mat % child("density")
@ -940,19 +802,10 @@ contains
ALL_NUCLIDES: do j = 1, mat % n_nuclides
! Check that this nuclide is listed in the cross_sections.xml file
name = trim(names % data(j))
if (.not. library_dict % has(to_lower(name))) then
if (.not. library_present(LIBRARY_NEUTRON, (to_lower(name)))) then
call fatal_error("Could not find nuclide " // trim(name) &
// " in cross_sections data file!")
end if
i_library = library_dict % get(to_lower(name))
if (run_CE) then
! Check to make sure cross-section is continuous energy neutron table
if (libraries(i_library) % type /= LIBRARY_NEUTRON) then
call fatal_error("Cross-section table " // trim(name) &
// " is not a continuous-energy neutron table.")
end if
end if
! If this nuclide hasn't been encountered yet, we need to add its name
! and alias to the nuclide_dict
@ -971,7 +824,7 @@ contains
element = name(1:scan(name, '0123456789') - 1)
! Make sure photon cross section data is available
if (.not. library_dict % has(to_lower(element))) then
if (.not. library_present(LIBRARY_PHOTON, to_lower(element))) then
call fatal_error("Could not find element " // trim(element) &
// " in cross_sections data file!")
end if
@ -1065,23 +918,11 @@ contains
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_dict % has(to_lower(name))) then
if (.not. library_present(LIBRARY_THERMAL, to_lower(name))) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
// " in cross_sections.xml file!")
end if
! Find index in xs_listing and set the name and alias according to the
! listing
i_library = library_dict % get(to_lower(name))
if (run_CE) then
! Check to make sure cross-section is continuous energy neutron table
if (libraries(i_library) % type /= LIBRARY_THERMAL) then
call fatal_error("Cross-section table " // trim(name) &
// " is not a S(a,b) table.")
end if
end if
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
if (.not. sab_dict % has(to_lower(name))) then
@ -2064,130 +1905,27 @@ contains
end subroutine read_plots_xml
!===============================================================================
! READ_*_CROSS_SECTIONS_XML reads information from a cross_sections.xml file. This
! file contains a listing of the CE and MG cross sections that may be used.
!===============================================================================
subroutine read_ce_cross_sections_xml()
subroutine read_mg_cross_sections_header() bind(C)
integer :: i ! loop index
integer :: n
integer :: n_libraries
logical :: file_exists ! does cross_sections.xml exist?
character(MAX_WORD_LEN) :: directory ! directory with cross sections
character(MAX_WORD_LEN) :: words(MAX_WORDS)
character(10000) :: temp_str
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_library
type(XMLNode), allocatable :: node_library_list(:)
! Check if cross_sections.xml exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Reading cross sections XML file...", 5)
! Parse cross_sections.xml file
call doc % load_file(path_cross_sections)
root = doc % document_element()
if (check_for_node(root, "directory")) then
! Copy directory information if present
call get_node_value(root, "directory", directory)
else
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
end if
! Get node list of all <library>
call get_node_list(root, "library", node_library_list)
n_libraries = size(node_library_list)
! Allocate xs_listings array
if (n_libraries == 0) then
call fatal_error("No cross section libraries present in cross_sections.xml &
&file!")
else
allocate(libraries(n_libraries))
end if
do i = 1, n_libraries
! Get pointer to ace table XML node
node_library = node_library_list(i)
! Get list of materials
if (check_for_node(node_library, "materials")) then
call get_node_value(node_library, "materials", temp_str)
call split_string(temp_str, words, n)
allocate(libraries(i) % materials(n))
libraries(i) % materials(:) = words(1:n)
end if
! Get type of library
if (check_for_node(node_library, "type")) then
call get_node_value(node_library, "type", temp_str)
select case(to_lower(temp_str))
case ('neutron')
libraries(i) % type = LIBRARY_NEUTRON
case ('thermal')
libraries(i) % type = LIBRARY_THERMAL
case ('photon')
libraries(i) % type = LIBRARY_PHOTON
end select
else
call fatal_error("Missing library type")
end if
! determine path of cross section table
if (check_for_node(node_library, "path")) then
call get_node_value(node_library, "path", temp_str)
else
call fatal_error("Missing library path")
end if
if (starts_with(temp_str, '/')) then
libraries(i) % path = trim(temp_str)
else
if (ends_with(directory,'/')) then
libraries(i) % path = trim(directory) // trim(temp_str)
else
libraries(i) % path = trim(directory) // '/' // trim(temp_str)
end if
end if
inquire(FILE=libraries(i) % path, EXIST=file_exists)
if (.not. file_exists) then
call warning("Cross section library " // trim(libraries(i) % path) // &
" does not exist.")
end if
end do
! Close cross sections XML file
call doc % clear()
end subroutine read_ce_cross_sections_xml
subroutine read_mg_cross_sections_header()
integer :: i ! loop index
integer :: n_libraries
logical :: file_exists ! does mgxs.h5 exist?
integer(HID_T) :: file_id
character(len=MAX_WORD_LEN), allocatable :: names(:)
character(kind=C_CHAR), pointer :: string(:)
interface
subroutine read_mg_cross_sections_header_c(file_id) bind(C)
import HID_T
integer(HID_T), value :: file_id
end subroutine
function path_cross_sections_c() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
end interface
call c_f_pointer(path_cross_sections_c(), string, [255])
path_cross_sections = to_f_string(string)
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
@ -2244,24 +1982,6 @@ contains
call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), &
energy_max(NEUTRON))
! Get the datasets present in the library
call get_groups(file_id, names)
n_libraries = size(names)
! Allocate libraries array
if (n_libraries == 0) then
call fatal_error("At least one MGXS data set must be present in &
&mgxs library file!")
else
allocate(libraries(n_libraries))
end if
do i = 1, n_libraries
! Get name of material
allocate(libraries(i) % materials(1))
libraries(i) % materials(1) = names(i)
end do
! Close MGXS HDF5 file
call file_close(file_id)
@ -2350,7 +2070,6 @@ contains
type(VectorReal), intent(in) :: sab_temps(:)
integer :: i, j
integer :: i_library
integer :: i_nuclide
integer :: i_element
integer :: i_sab
@ -2358,6 +2077,7 @@ contains
integer(HID_T) :: group_id
logical :: mp_found ! if windowed multipole libraries were found
character(MAX_WORD_LEN) :: name
character(MAX_FILE_LEN) :: filename
character(3) :: element
type(SetChar) :: already_read
type(SetChar) :: element_already_read
@ -2375,14 +2095,14 @@ contains
name = materials(i) % names(j)
if (.not. already_read % contains(name)) then
i_library = library_dict % get(to_lower(name))
filename = library_path(LIBRARY_NEUTRON, to_lower(name))
i_nuclide = nuclide_dict % get(to_lower(name))
call write_message('Reading ' // trim(name) // ' from ' // &
trim(libraries(i_library) % path), 6)
trim(filename), 6)
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read nuclide data from HDF5
@ -2415,13 +2135,13 @@ contains
if (photon_transport) then
if (.not. element_already_read % contains(element)) then
! Read photon interaction data from HDF5 photon library
i_library = library_dict % get(to_lower(element))
filename = library_path(LIBRARY_PHOTON, to_lower(element))
i_element = element_dict % get(element)
call write_message('Reading ' // trim(element) // ' from ' // &
trim(libraries(i_library) % path), 6)
trim(filename), 6)
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read element data from HDF5
@ -2505,14 +2225,14 @@ contains
name = materials(i) % sab_names(j)
if (.not. already_read % contains(name)) then
i_library = library_dict % get(to_lower(name))
filename = library_path(LIBRARY_THERMAL, to_lower(name))
i_sab = sab_dict % get(to_lower(name))
call write_message('Reading ' // trim(name) // ' from ' // &
trim(libraries(i_library) % path), 6)
trim(filename), 6)
! Open file and make sure version matches
file_id = file_open(libraries(i_library) % path, 'r')
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read S(a,b) data from HDF5
@ -2556,9 +2276,8 @@ contains
end if
end do
if (.not. mp_found) call warning("Windowed multipole functionality is &
&turned on, but no multipole libraries were found. Set the &
&<multipole_library> element in settings.xml or the &
&OPENMC_MULTIPOLE_LIBRARY environment variable.")
&turned on, but no multipole libraries were found. Make sure that &
&windowed multipole data is present in your cross_sections.xml file.")
end if
call already_read % clear()
@ -2578,20 +2297,18 @@ contains
logical :: file_exists ! Does multipole library exist?
character(7) :: readable ! Is multipole library readable?
character(MAX_FILE_LEN) :: filename ! Path to multipole xs library
character(kind=C_CHAR), pointer :: string(:)
integer(HID_T) :: file_id
integer(HID_T) :: group_id
! For the time being, and I know this is a bit hacky, we just assume
! that the file will be ZZZAAAmM.h5.
associate (nuc => nuclides(i_table))
if (nuc % metastable > 0) then
filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // &
trim(zero_padded(nuc % A, 3)) // 'm' // &
trim(to_str(nuc % metastable)) // ".h5"
! Look for WMP data in cross_sections.xml
if (library_present(LIBRARY_WMP, to_lower(nuc % name))) then
filename = library_path(LIBRARY_WMP, to_lower(nuc % name))
else
filename = trim(path_multipole) // trim(zero_padded(nuc % Z, 3)) // &
trim(zero_padded(nuc % A, 3)) // ".h5"
nuc % mp_present = .false.
return
end if
! Check if Multipole library exists and is readable
@ -2605,7 +2322,8 @@ contains
end if
! Display message
call write_message("Loading Windowed Multipole XS from " // filename, 6)
call write_message("Reading " // trim(nuc % name) // " WMP data from " &
// filename, 6)
! Open file and make sure version is sufficient
file_id = file_open(filename, 'r')

37
src/lattice.cpp
View file

@ -18,10 +18,13 @@ namespace openmc {
// Global variables
//==============================================================================
std::vector<Lattice*> lattices;
namespace model {
std::vector<Lattice*> lattices;
std::unordered_map<int32_t, int32_t> lattice_map;
}
//==============================================================================
// Lattice implementation
//==============================================================================
@ -65,8 +68,8 @@ Lattice::adjust_indices()
// Adjust the indices for the universes array.
for (LatticeIter it = begin(); it != end(); ++it) {
int uid = *it;
auto search = universe_map.find(uid);
if (search != universe_map.end()) {
auto search = model::universe_map.find(uid);
if (search != model::universe_map.end()) {
*it = search->second;
} else {
std::stringstream err_msg;
@ -78,8 +81,8 @@ Lattice::adjust_indices()
// Adjust the index for the outer universe.
if (outer_ != NO_OUTER_UNIVERSE) {
auto search = universe_map.find(outer_);
if (search != universe_map.end()) {
auto search = model::universe_map.find(outer_);
if (search != model::universe_map.end()) {
outer_ = search->second;
} else {
std::stringstream err_msg;
@ -118,7 +121,7 @@ Lattice::to_hdf5(hid_t lattices_group) const
}
if (outer_ != NO_OUTER_UNIVERSE) {
int32_t outer_id = universes[outer_]->id_;
int32_t outer_id = model::universes[outer_]->id_;
write_dataset(lat_group, "outer", outer_id);
} else {
write_dataset(lat_group, "outer", outer_);
@ -370,7 +373,7 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const
for (int j = 0; j < nx; j++) {
int indx1 = nx*ny*m + nx*k + j;
int indx2 = nx*ny*m + nx*(ny-k-1) + j;
out[indx2] = universes[universes_[indx1]]->id_;
out[indx2] = model::universes[universes_[indx1]]->id_;
}
}
}
@ -387,7 +390,7 @@ RectLattice::to_hdf5_inner(hid_t lat_group) const
for (int j = 0; j < nx; j++) {
int indx1 = nx*k + j;
int indx2 = nx*(ny-k-1) + j;
out[indx2] = universes[universes_[indx1]]->id_;
out[indx2] = model::universes[universes_[indx1]]->id_;
}
}
@ -847,7 +850,7 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
// This array position is never used; put a -1 to indicate this.
out[indx] = -1;
} else {
out[indx] = universes[universes_[indx]]->id_;
out[indx] = model::universes[universes_[indx]]->id_;
}
}
}
@ -865,18 +868,18 @@ extern "C" void
read_lattices(pugi::xml_node *node)
{
for (pugi::xml_node lat_node : node->children("lattice")) {
lattices.push_back(new RectLattice(lat_node));
model::lattices.push_back(new RectLattice(lat_node));
}
for (pugi::xml_node lat_node : node->children("hex_lattice")) {
lattices.push_back(new HexLattice(lat_node));
model::lattices.push_back(new HexLattice(lat_node));
}
// Fill the lattice map.
for (int i_lat = 0; i_lat < lattices.size(); i_lat++) {
int id = lattices[i_lat]->id_;
auto in_map = lattice_map.find(id);
if (in_map == lattice_map.end()) {
lattice_map[id] = i_lat;
for (int i_lat = 0; i_lat < model::lattices.size(); i_lat++) {
int id = model::lattices[i_lat]->id_;
auto in_map = model::lattice_map.find(id);
if (in_map == model::lattice_map.end()) {
model::lattice_map[id] = i_lat;
} else {
std::stringstream err_msg;
err_msg << "Two or more lattices use the same unique ID: " << id;
@ -890,7 +893,7 @@ read_lattices(pugi::xml_node *node)
//==============================================================================
extern "C" {
Lattice* lattice_pointer(int lat_ind) {return lattices[lat_ind];}
Lattice* lattice_pointer(int lat_ind) {return model::lattices[lat_ind];}
int32_t lattice_id(Lattice *lat) {return lat->id_;}
}

24
src/main.cpp
View file

@ -1,13 +1,15 @@
#ifdef OPENMC_MPI
#include "mpi.h"
#include <mpi.h>
#endif
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
int main(int argc, char* argv[]) {
using namespace openmc;
int err;
// Initialize run -- when run with MPI, pass communicator
@ -21,30 +23,30 @@ int main(int argc, char* argv[]) {
// This happens for the -h and -v flags
return 0;
} else if (err) {
openmc::fatal_error(openmc_err_msg);
fatal_error(openmc_err_msg);
}
// start problem based on mode
switch (openmc::settings::run_mode) {
case openmc::RUN_MODE_FIXEDSOURCE:
case openmc::RUN_MODE_EIGENVALUE:
switch (settings::run_mode) {
case RUN_MODE_FIXEDSOURCE:
case RUN_MODE_EIGENVALUE:
err = openmc_run();
break;
case openmc::RUN_MODE_PLOTTING:
case RUN_MODE_PLOTTING:
err = openmc_plot_geometry();
break;
case openmc::RUN_MODE_PARTICLE:
if (openmc_master) err = openmc_particle_restart();
case RUN_MODE_PARTICLE:
if (mpi::master) err = openmc_particle_restart();
break;
case openmc::RUN_MODE_VOLUME:
case RUN_MODE_VOLUME:
err = openmc_calculate_volumes();
break;
}
if (err) openmc::fatal_error(openmc_err_msg);
if (err) fatal_error(openmc_err_msg);
// Finalize and free up memory
err = openmc_finalize();
if (err) openmc::fatal_error(openmc_err_msg);
if (err) fatal_error(openmc_err_msg);
// If MPI is in use and enabled, terminate it
#ifdef OPENMC_MPI

40
src/material.cpp
View file

@ -13,9 +13,13 @@ namespace openmc {
// Global variables
//==============================================================================
namespace model {
std::vector<Material*> materials;
std::unordered_map<int32_t, int32_t> material_map;
} // namespace model
//==============================================================================
// Material implementation
//==============================================================================
@ -46,16 +50,16 @@ read_materials(pugi::xml_node* node)
{
// Loop over XML material elements and populate the array.
for (pugi::xml_node material_node : node->children("material")) {
materials.push_back(new Material(material_node));
model::materials.push_back(new Material(material_node));
}
materials.shrink_to_fit();
model::materials.shrink_to_fit();
// Populate the material map.
for (int i = 0; i < materials.size(); i++) {
int32_t mid = materials[i]->id_;
auto search = material_map.find(mid);
if (search == material_map.end()) {
material_map[mid] = i;
for (int i = 0; i < model::materials.size(); i++) {
int32_t mid = model::materials[i]->id_;
auto search = model::material_map.find(mid);
if (search == model::material_map.end()) {
model::material_map[mid] = i;
} else {
std::stringstream err_msg;
err_msg << "Two or more materials use the same unique ID: " << mid;
@ -71,8 +75,8 @@ read_materials(pugi::xml_node* node)
extern "C" int
openmc_material_get_volume(int32_t index, double* volume)
{
if (index >= 1 && index <= materials.size()) {
Material* m = materials[index - 1];
if (index >= 1 && index <= model::materials.size()) {
Material* m = model::materials[index - 1];
if (m->volume_ >= 0.0) {
*volume = m->volume_;
return 0;
@ -91,8 +95,8 @@ openmc_material_get_volume(int32_t index, double* volume)
extern "C" int
openmc_material_set_volume(int32_t index, double volume)
{
if (index >= 1 && index <= materials.size()) {
Material* m = materials[index - 1];
if (index >= 1 && index <= model::materials.size()) {
Material* m = model::materials[index - 1];
if (volume >= 0.0) {
m->volume_ = volume;
return 0;
@ -111,7 +115,7 @@ openmc_material_set_volume(int32_t index, double volume)
//==============================================================================
extern "C" {
Material* material_pointer(int32_t indx) {return materials[indx];}
Material* material_pointer(int32_t indx) {return model::materials[indx];}
int32_t material_id(Material* mat) {return mat->id_;}
@ -119,7 +123,7 @@ extern "C" {
{
mat->id_ = id;
//TODO: off-by-one
material_map[id] = index - 1;
model::material_map[id] = index - 1;
}
bool material_fissionable(Material* mat) {return mat->fissionable;}
@ -131,17 +135,17 @@ extern "C" {
void extend_materials_c(int32_t n)
{
materials.reserve(materials.size() + n);
model::materials.reserve(model::materials.size() + n);
for (int32_t i = 0; i < n; i++) {
materials.push_back(new Material());
model::materials.push_back(new Material());
}
}
void free_memory_material_c()
{
for (Material *mat : materials) {delete mat;}
materials.clear();
material_map.clear();
for (Material *mat : model::materials) {delete mat;}
model::materials.clear();
model::material_map.clear();
}
}

70
src/material_header.F90
View file

@ -12,7 +12,7 @@ module material_header
use sab_header
use simulation_header, only: log_spacing
use stl_vector, only: VectorReal, VectorInt
use string, only: to_str
use string, only: to_str, to_f_string
implicit none
@ -166,34 +166,52 @@ contains
call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL))
end subroutine material_set_fissionable
function material_set_density(this, density) result(err)
function material_set_density(this, density, units) result(err)
class(Material), intent(inout) :: this
real(8), intent(in) :: density
character(*), intent(in) :: units
integer :: err
integer :: i
real(8) :: sum_percent
real(8) :: awr
real(8) :: previous_density_gpcc
real(8) :: f
if (allocated(this % atom_density)) then
! Set total density based on value provided
this % density = density
! Determine normalized atom percents
sum_percent = sum(this % atom_density)
this % atom_density(:) = this % atom_density / sum_percent
! Recalculate nuclide atom densities based on given density
this % atom_density(:) = density * this % atom_density
! Calculate density in g/cm^3.
this % density_gpcc = ZERO
do i = 1, this % n_nuclides
awr = nuclides(this % nuclide(i)) % awr
this % density_gpcc = this % density_gpcc &
+ this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
end do
err = 0
select case (units)
case ('atom/b-cm')
! Set total density based on value provided
this % density = density
! Determine normalized atom percents
sum_percent = sum(this % atom_density)
this % atom_density(:) = this % atom_density / sum_percent
! Recalculate nuclide atom densities based on given density
this % atom_density(:) = density * this % atom_density
! Calculate density in g/cm^3.
this % density_gpcc = ZERO
do i = 1, this % n_nuclides
awr = nuclides(this % nuclide(i)) % awr
this % density_gpcc = this % density_gpcc &
+ this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
end do
case ('g/cm3', 'g/cc')
! Determine factor by which to change densities
previous_density_gpcc = this % density_gpcc
f = density / previous_density_gpcc
! Update densities
this % density_gpcc = density
this % density = f * this % density
this % atom_density(:) = f * this % atom_density(:)
case default
err = E_INVALID_ARGUMENT
call set_errmsg("Invalid units '" // trim(units) // "' specified.")
end select
else
err = E_ALLOCATE
call set_errmsg("Material atom density array hasn't been allocated.")
@ -738,16 +756,22 @@ contains
end function openmc_material_set_id
function openmc_material_set_density(index, density) result(err) bind(C)
! Set the total density of a material in atom/b-cm
function openmc_material_set_density(index, density, units) result(err) bind(C)
! Set the total density of a material
integer(C_INT32_T), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: density
character(kind=C_CHAR), intent(in) :: units(*)
integer(C_INT) :: err
character(:), allocatable :: units_
! Convert C string to Fortran string
units_ = to_f_string(units)
err = E_UNASSIGNED
if (index >= 1 .and. index <= size(materials)) then
associate (m => materials(index))
err = m % set_density(density)
err = m % set_density(density, units_)
end associate
else
err = E_OUT_OF_BOUNDS
@ -794,7 +818,7 @@ contains
m % n_nuclides = n
! Set total density to the sum of the vector
err = m % set_density(sum(density))
err = m % set_density(sum(density), 'atom/b-cm')
! Assign S(a,b) tables
call m % assign_sab_tables()

61
src/mesh.cpp
View file

@ -29,10 +29,13 @@ namespace openmc {
// Global variables
//==============================================================================
std::vector<std::unique_ptr<RegularMesh>> meshes;
namespace model {
std::vector<std::unique_ptr<RegularMesh>> meshes;
std::unordered_map<int32_t, int32_t> mesh_map;
} // namespace model
//==============================================================================
// RegularMesh implementation
//==============================================================================
@ -44,7 +47,7 @@ RegularMesh::RegularMesh(pugi::xml_node node)
id_ = std::stoi(get_node_value(node, "id"));
// Check to make sure 'id' hasn't been used
if (mesh_map.find(id_) != mesh_map.end()) {
if (model::mesh_map.find(id_) != model::mesh_map.end()) {
fatal_error("Two or more meshes use the same unique ID: " +
std::to_string(id_));
}
@ -724,11 +727,11 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
extern "C" int
openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
{
if (index_start) *index_start = meshes.size();
if (index_start) *index_start = model::meshes.size();
for (int i = 0; i < n; ++i) {
meshes.emplace_back(new RegularMesh{});
model::meshes.emplace_back(new RegularMesh{});
}
if (index_end) *index_end = meshes.size() - 1;
if (index_end) *index_end = model::meshes.size() - 1;
return 0;
}
@ -737,8 +740,8 @@ openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
extern "C" int
openmc_get_mesh_index(int32_t id, int32_t* index)
{
auto pair = mesh_map.find(id);
if (pair == mesh_map.end()) {
auto pair = model::mesh_map.find(id);
if (pair == model::mesh_map.end()) {
set_errmsg("No mesh exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
}
@ -750,11 +753,11 @@ openmc_get_mesh_index(int32_t id, int32_t* index)
extern "C" int
openmc_mesh_get_id(int32_t index, int32_t* id)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*id = meshes[index]->id_;
*id = model::meshes[index]->id_;
return 0;
}
@ -762,12 +765,12 @@ openmc_mesh_get_id(int32_t index, int32_t* id)
extern "C" int
openmc_mesh_set_id(int32_t index, int32_t id)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
meshes[index]->id_ = id;
mesh_map[id] = index;
model::meshes[index]->id_ = id;
model::mesh_map[id] = index;
return 0;
}
@ -775,12 +778,12 @@ openmc_mesh_set_id(int32_t index, int32_t id)
extern "C" int
openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
*dims = meshes[index]->shape_.data();
*n = meshes[index]->n_dimension_;
*dims = model::meshes[index]->shape_.data();
*n = model::meshes[index]->n_dimension_;
return 0;
}
@ -788,14 +791,14 @@ openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
extern "C" int
openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// Copy dimension
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
auto& m = meshes[index];
auto& m = model::meshes[index];
m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
m->n_dimension_ = m->shape_.size();
@ -806,12 +809,12 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
extern "C" int
openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
auto& m = meshes[index];
auto& m = model::meshes[index];
if (m->lower_left_.dimension() == 0) {
set_errmsg("Mesh parameters have not been set.");
return OPENMC_E_ALLOCATE;
@ -829,12 +832,12 @@ extern "C" int
openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
const double* width)
{
if (index < 0 || index >= meshes.size()) {
if (index < 0 || index >= model::meshes.size()) {
set_errmsg("Index in meshes array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
auto& m = meshes[index];
auto& m = model::meshes[index];
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
if (ll && ur) {
m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
@ -864,10 +867,10 @@ void read_meshes(pugi::xml_node* root)
{
for (auto node : root->children("mesh")) {
// Read mesh and add to vector
meshes.emplace_back(new RegularMesh{node});
model::meshes.emplace_back(new RegularMesh{node});
// Map ID to position in vector
mesh_map[meshes.back()->id_] = meshes.size() - 1;
model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
}
}
@ -875,13 +878,13 @@ void meshes_to_hdf5(hid_t group)
{
// Write number of meshes
hid_t meshes_group = create_group(group, "meshes");
int32_t n_meshes = meshes.size();
int32_t n_meshes = model::meshes.size();
write_attribute(meshes_group, "n_meshes", n_meshes);
if (n_meshes > 0) {
// Write IDs of meshes
std::vector<int> ids;
for (const auto& m : meshes) {
for (const auto& m : model::meshes) {
m->to_hdf5(meshes_group);
ids.push_back(m->id_);
}
@ -896,9 +899,9 @@ void meshes_to_hdf5(hid_t group)
//==============================================================================
extern "C" {
int n_meshes() { return meshes.size(); }
int n_meshes() { return model::meshes.size(); }
RegularMesh* mesh_ptr(int i) { return meshes.at(i).get(); }
RegularMesh* mesh_ptr(int i) { return model::meshes.at(i).get(); }
int32_t mesh_id(RegularMesh* m) { return m->id_; }
@ -936,8 +939,8 @@ extern "C" {
void free_memory_mesh()
{
meshes.clear();
mesh_map.clear();
model::meshes.clear();
model::mesh_map.clear();
}
}

10
src/mgxs.cpp
View file

@ -17,15 +17,23 @@
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/string_functions.h"
#include "openmc/string_utils.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace data {
// Storage for the MGXS data
std::vector<Mgxs> nuclides_MG;
std::vector<Mgxs> macro_xs;
} // namespace data
//==============================================================================
// Mgxs base-class methods
//==============================================================================

57
src/mgxs_interface.cpp
View file

@ -2,6 +2,7 @@
#include <string>
#include "openmc/cross_sections.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -12,10 +13,14 @@ namespace openmc {
// Global variable definitions
//==============================================================================
namespace data {
std::vector<double> energy_bins;
std::vector<double> energy_bin_avg;
std::vector<double> rev_energy_bins;
} // namesapce data
//==============================================================================
// Mgxs data loading interface methods
//==============================================================================
@ -43,7 +48,7 @@ add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance,
max_order, legendre_to_tabular, legendre_to_tabular_points, method);
nuclides_MG.push_back(mg);
data::nuclides_MG.push_back(mg);
close_group(xs_grp);
}
@ -54,7 +59,7 @@ query_fissionable_c(int n_nuclides, const int i_nuclides[])
{
bool result = false;
for (int n = 0; n < n_nuclides; n++) {
if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
if (data::nuclides_MG[i_nuclides[n] - 1].fissionable) result = true;
}
return result;
}
@ -78,16 +83,16 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
// material
std::vector<Mgxs*> mgxs_ptr(n_nuclides);
for (int n = 0; n < n_nuclides; n++) {
mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
mgxs_ptr[n] = &data::nuclides_MG[i_nuclides[n] - 1];
}
Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
tolerance, method);
macro_xs.emplace_back(macro);
data::macro_xs.emplace_back(macro);
} else {
// Preserve the ordering of materials by including a blank entry
Mgxs macro;
macro_xs.emplace_back(macro);
data::macro_xs.emplace_back(macro);
}
}
@ -96,18 +101,32 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
void read_mg_cross_sections_header_c(hid_t file_id)
{
ensure_exists(file_id, "energy_groups", true);
read_attribute(file_id, "energy_groups", num_energy_groups);
read_attribute(file_id, "energy_groups", data::num_energy_groups);
ensure_exists(file_id, "group structure", true);
read_attribute(file_id, "group structure", rev_energy_bins);
read_attribute(file_id, "group structure", data::rev_energy_bins);
// Reverse energy bins
std::copy(rev_energy_bins.crbegin(), rev_energy_bins.crend(),
std::back_inserter(energy_bins));
std::copy(data::rev_energy_bins.crbegin(), data::rev_energy_bins.crend(),
std::back_inserter(data::energy_bins));
// Create average energies
for (int i = 0; i < energy_bins.size() - 1; ++i) {
energy_bin_avg.push_back(0.5*(energy_bins[i] + energy_bins[i+1]));
for (int i = 0; i < data::energy_bins.size() - 1; ++i) {
data::energy_bin_avg.push_back(0.5*(data::energy_bins[i] + data::energy_bins[i+1]));
}
// Add entries into libraries for MG data
auto names = group_names(file_id);
if (names.empty()) {
fatal_error("At least one MGXS data set must be present in mgxs "
"library file!");
}
for (auto& name : names) {
Library lib {};
lib.type_ = Library::Type::neutron;
lib.materials_.push_back(name);
data::libraries.push_back(lib);
}
}
@ -119,7 +138,7 @@ void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
data::macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
nu_fiss_xs);
}
@ -144,7 +163,7 @@ get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
} else {
dg_c_p = dg;
}
return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
return data::nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
@ -168,7 +187,7 @@ get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
} else {
dg_c_p = dg;
}
return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
@ -177,7 +196,7 @@ void
set_nuclide_angle_index_c(int index, const double uvw[3])
{
// Update the values
nuclides_MG[index - 1].set_angle_index(uvw);
data::nuclides_MG[index - 1].set_angle_index(uvw);
}
//==============================================================================
@ -186,7 +205,7 @@ void
set_macro_angle_index_c(int index, const double uvw[3])
{
// Update the values
macro_xs[index - 1].set_angle_index(uvw);
data::macro_xs[index - 1].set_angle_index(uvw);
}
//==============================================================================
@ -195,7 +214,7 @@ void
set_nuclide_temperature_index_c(int index, double sqrtkT)
{
// Update the values
nuclides_MG[index - 1].set_temperature_index(sqrtkT);
data::nuclides_MG[index - 1].set_temperature_index(sqrtkT);
}
//==============================================================================
@ -210,7 +229,7 @@ get_name_c(int index, int name_len, char* name)
std::strcpy(name, str.c_str());
// Now get the data and copy to the C-string
str = nuclides_MG[index - 1].name;
str = data::nuclides_MG[index - 1].name;
std::strcpy(name, str.c_str());
// Finally, remove the null terminator
@ -222,7 +241,7 @@ get_name_c(int index, int name_len, char* name)
double
get_awr_c(int index)
{
return nuclides_MG[index - 1].awr;
return data::nuclides_MG[index - 1].awr;
}
} // namespace openmc

8
src/nuclide.cpp
View file

@ -6,9 +6,13 @@ namespace openmc {
// Global variables
//==============================================================================
namespace data {
std::array<double, 2> energy_min {0.0, 0.0};
std::array<double, 2> energy_max {INFTY, INFTY};
} // namespace data
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -16,8 +20,8 @@ std::array<double, 2> energy_max {INFTY, INFTY};
extern "C" void
set_particle_energy_bounds(int particle, double E_min, double E_max)
{
energy_min[particle - 1] = E_min;
energy_max[particle - 1] = E_max;
data::energy_min[particle - 1] = E_min;
data::energy_max[particle - 1] = E_max;
}
} // namespace openmc

118
src/nuclide_header.F90
View file

@ -6,7 +6,8 @@ module nuclide_header
use algorithm, only: sort, find, binary_search
use constants
use dict_header, only: DictIntInt, DictCharInt
use endf, only: reaction_name, is_fission, is_disappearance
use endf, only: reaction_name, is_fission, is_disappearance, &
is_inelastic_scatter
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error
use hdf5_interface
@ -176,20 +177,6 @@ module nuclide_header
real(C_DOUBLE) :: pair_production ! macroscopic pair production xs
end type MaterialMacroXS
!===============================================================================
! LIBRARY contains data read from a cross_sections.xml file
!===============================================================================
type Library
integer :: type
character(MAX_WORD_LEN), allocatable :: materials(:)
character(MAX_FILE_LEN) :: path
end type Library
! Cross section libraries
type(Library), allocatable :: libraries(:)
type(DictCharInt) :: library_dict
! Nuclear data for each nuclide
type(Nuclide), allocatable, target :: nuclides(:)
integer(C_INT), bind(C) :: n_nuclides
@ -204,8 +191,46 @@ module nuclide_header
real(8) :: energy_min(2) = [ZERO, ZERO]
real(8) :: energy_max(2) = [INFINITY, INFINITY]
interface
function library_present_c(type, name) result(b) bind(C, name='library_present')
import C_INT, C_CHAR, C_BOOL
integer(C_INT), value :: type
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: b
end function
function library_path_c(type, name) result(path) bind(C, name='library_path')
import C_INT, C_CHAR, C_PTR
integer(C_INT), value :: type
character(kind=C_CHAR), intent(in) :: name(*)
type(C_PTR) :: path
end function
end interface
contains
function library_path(type, name) result(path)
integer, intent(in) :: type
character(len=*), intent(in) :: name
character(MAX_FILE_LEN) :: path
type(C_PTR) :: ptr
character(kind=C_CHAR), pointer :: string(:)
ptr = library_path_c(type, to_c_string(name))
call c_f_pointer(ptr, string, [255])
path = to_f_string(string)
end function
function library_present(type, name) result(b)
integer, intent(in) :: type
character(len=*), intent(in) :: name
logical :: b
b = library_present_c(type, to_c_string(name))
end function
!===============================================================================
! ASSIGN_0K_ELASTIC_SCATTERING
!===============================================================================
@ -216,11 +241,30 @@ contains
integer :: i
real(8) :: xs_cdf_sum
interface
function res_scat_nuclides_empty() result(empty) bind(C)
import C_BOOL
logical(C_BOOL) :: empty
end function
function res_scat_nuclides_size() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
function res_scat_nuclides_cmp(i, name) result(b) bind(C)
import C_INT, C_CHAR, C_BOOL
integer(C_INT), value :: i
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: b
end function
end interface
this % resonant = .false.
if (allocated(res_scat_nuclides)) then
if (.not. res_scat_nuclides_empty()) then
! If resonant nuclides were specified, check the list explicitly
do i = 1, size(res_scat_nuclides)
if (this % name == res_scat_nuclides(i)) then
do i = 1, res_scat_nuclides_size()
if (res_scat_nuclides_cmp(i, to_c_string(this % name))) then
this % resonant = .true.
! Make sure nuclide has 0K data
@ -486,12 +530,7 @@ contains
! Add the reaction index to the scattering array if this is an inelastic
! scatter reaction
if (MTs % data(i) /= N_FISSION .and. MTs % data(i) /= N_F .and. &
MTs % data(i) /= N_NF .and. MTs % data(i) /= N_2NF .and. &
MTs % data(i) /= N_3NF .and. MTs % data(i) < 200 .and. &
MTs % data(i) /= N_LEVEL .and. MTs % data(i) /= ELASTIC .and. &
.not. this % reactions(i) % redundant) then
if (is_inelastic_scatter(MTs % data(i))) then
call index_inelastic_scatter % push_back(i)
end if
@ -653,20 +692,10 @@ contains
end if
end do
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rx % MT == N_LEVEL .or. rx % MT == N_NONELASTIC) cycle
! Skip gas production cross sections (MT=200+), etc.
if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle
! Skip level cross sections if total is available
if (rx % MT >= N_P0 .and. rx % MT <= N_PC .and. find(MTs, N_P) /= -1) cycle
if (rx % MT >= N_D0 .and. rx % MT <= N_DC .and. find(MTs, N_D) /= -1) cycle
if (rx % MT >= N_T0 .and. rx % MT <= N_TC .and. find(MTs, N_T) /= -1) cycle
if (rx % MT >= N_3HE0 .and. rx % MT <= N_3HEC .and. find(MTs, N_3HE) /= -1) cycle
if (rx % MT >= N_A0 .and. rx % MT <= N_AC .and. find(MTs, N_A) /= -1) cycle
if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle
! Skip redundant reactions, which are used for photon production
! Skip any reaction that has been marked as redundant
if (rx % redundant) cycle
! Add contribution to total cross section
@ -1540,6 +1569,11 @@ contains
subroutine free_memory_nuclide()
integer :: i
interface
subroutine library_clear() bind(C)
end subroutine
end interface
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) then
! First call the clear routines
@ -1550,10 +1584,8 @@ contains
end if
n_nuclides = 0
if (allocated(libraries)) deallocate(libraries)
call nuclide_dict % clear()
call library_dict % clear()
call library_clear()
end subroutine free_memory_nuclide
@ -1593,11 +1625,11 @@ contains
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT) :: err
integer :: i_library
integer :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
character(:), allocatable :: name_
character(MAX_FILE_LEN) :: filename
real(8) :: minmax(2) = [ZERO, INFINITY]
type(VectorReal) :: temperature
type(Nuclide), allocatable :: new_nuclides(:)
@ -1607,7 +1639,7 @@ contains
err = 0
if (.not. nuclide_dict % has(to_lower(name_))) then
if (library_dict % has(to_lower(name_))) then
if (library_present(LIBRARY_NEUTRON, to_lower(name_))) then
! allocate extra space in nuclides array
n = n_nuclides
allocate(new_nuclides(n + 1))
@ -1615,10 +1647,10 @@ contains
call move_alloc(FROM=new_nuclides, TO=nuclides)
n = n + 1
i_library = library_dict % get(to_lower(name_))
filename = library_path(LIBRARY_NEUTRON, to_lower(name_))
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read nuclide data from HDF5

32
src/output.cpp
View file

@ -51,7 +51,7 @@ std::string time_stamp()
<< ":" << now->tm_min << ":" << now->tm_sec;
return ts.str();
}
//==============================================================================
//===============================================================================
@ -62,7 +62,7 @@ void print_plot() {
header("PLOTTING SUMMARY", 5);
for (auto pl : plots) {
for (auto pl : model::plots) {
// Plot id
std::cout << "Plot ID: " << pl.id_ << "\n";
// Plot filename
@ -74,7 +74,7 @@ void print_plot() {
if (PlotType::slice == pl.type_) {
std::cout << "Plot Type: Slice" << "\n";
} else if (PlotType::voxel == pl.type_) {
std::cout << "Plot Type: Voxel" << "\n";
std::cout << "Plot Type: Voxel" << "\n";
}
// Plot parameters
@ -98,9 +98,9 @@ void print_plot() {
if (PlotColorBy::cells == pl.color_by_) {
std::cout << "Coloring: Cells" << "\n";
} else if (PlotColorBy::mats == pl.color_by_) {
std::cout << "Coloring: Materials" << "\n";
std::cout << "Coloring: Materials" << "\n";
}
if (PlotType::slice == pl.type_) {
switch(pl.basis_) {
case PlotBasis::xy:
@ -122,29 +122,29 @@ void print_plot() {
}
std::cout << "\n";
}
}
void
print_overlap_check() {
#ifdef OPENMC_MPI
std::vector<int64_t> temp(overlap_check_count);
int err = MPI_Reduce(temp.data(), overlap_check_count.data(),
overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0,
std::vector<int64_t> temp(model::overlap_check_count);
int err = MPI_Reduce(temp.data(), model::overlap_check_count.data(),
model::overlap_check_count.size(), MPI_INT64_T, MPI_SUM, 0,
mpi::intracomm);
#endif
if (openmc_master) {
if (mpi::master) {
header("cell overlap check summary", 1);
std::cout << " Cell ID No. Overlap Checks\n";
std::vector<int32_t> sparse_cell_ids;
for (int i = 0; i < n_cells; i++) {
std::cout << " " << std::setw(8) << cells[i]->id_ << std::setw(17)
<< overlap_check_count[i] << "\n";
if (overlap_check_count[i] < 10) {
sparse_cell_ids.push_back(cells[i]->id_);
for (int i = 0; i < model::cells.size(); i++) {
std::cout << " " << std::setw(8) << model::cells[i]->id_ << std::setw(17)
<< model::overlap_check_count[i] << "\n";
if (model::overlap_check_count[i] < 10) {
sparse_cell_ids.push_back(model::cells[i]->id_);
}
}

2
src/particle.cpp
View file

@ -113,7 +113,7 @@ Particle::from_source(const Bank* src)
} else {
g = static_cast<int>(src->E);
last_g = static_cast<int>(src->E);
E = energy_bin_avg[g - 1];
E = data::energy_bin_avg[g - 1];
}
last_E = E;
}

22
src/physics_mg.cpp
View file

@ -21,14 +21,13 @@
namespace openmc {
void
collision_mg(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs)
collision_mg(Particle* p, const MaterialMacroXS* material_xs)
{
// Add to the collision counter for the particle
p->n_collision++;
// Sample the reaction type
sample_reaction(p, energy_bin_avg, material_xs);
sample_reaction(p, material_xs);
// Display information about collision
if ((settings::verbosity >= 10) || (simulation::trace)) {
@ -39,21 +38,20 @@ collision_mg(Particle* p, const double* energy_bin_avg,
}
void
sample_reaction(Particle* p, const double* energy_bin_avg,
const MaterialMacroXS* material_xs)
sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
{
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
// change when sampling fission sites. The following block handles all
// absorption (including fission)
if (materials[p->material - 1]->fissionable) {
if (model::materials[p->material - 1]->fissionable) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
Bank* result_bank;
int64_t result_bank_size;
// Get pointer to fission bank from Fortran side
openmc_fission_bank(&result_bank, &result_bank_size);
create_fission_sites(p, result_bank, &n_bank, result_bank_size,
create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size,
material_xs);
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
@ -72,7 +70,7 @@ sample_reaction(Particle* p, const double* energy_bin_avg,
if (!p->alive) return;
// Sample a scattering event to determine the energy of the exiting neutron
scatter(p, energy_bin_avg);
scatter(p);
// Play Russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
@ -82,14 +80,14 @@ sample_reaction(Particle* p, const double* energy_bin_avg,
}
void
scatter(Particle* p, const double* energy_bin_avg)
scatter(Particle* p)
{
// Adjust indices for Fortran to C++ indexing
// TODO: Remove when no longer needed
int gin = p->last_g - 1;
int gout = p->g - 1;
int i_mat = p->material - 1;
macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
// Adjust return value for fortran indexing
// TODO: Remove when no longer needed
@ -99,7 +97,7 @@ scatter(Particle* p, const double* energy_bin_avg)
rotate_angle_c(p->coord[0].uvw, p->mu, nullptr);
// Update energy value for downstream compatability (in tallying)
p->E = energy_bin_avg[gout];
p->E = data::energy_bin_avg[gout];
// Set event component
p->event = EVENT_SCATTER;
@ -180,7 +178,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
// the energy in the fission bank
int dg;
int gout;
macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
data::macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
bank_array[i].E = static_cast<double>(gout + 1);
bank_array[i].delayed_group = dg + 1;

77
src/plot.cpp
View file

@ -9,7 +9,8 @@
#include "openmc/geometry.h"
#include "openmc/cell.h"
#include "openmc/material.h"
#include "openmc/string_functions.h"
#include "openmc/message_passing.h"
#include "openmc/string_utils.h"
#include "openmc/mesh.h"
#include "openmc/output.h"
#include "openmc/hdf5_interface.h"
@ -19,19 +20,23 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
const RGBColor WHITE {255, 255, 255};
constexpr int PLOT_LEVEL_LOWEST {-1}; //!< lower bound on plot universe level
//==============================================================================
// Global variables
//==============================================================================
int PLOT_LEVEL_LOWEST = -1;
std::unordered_map<int, int> plot_map;
int n_plots;
namespace model {
std::vector<Plot> plots;
std::unordered_map<int, int> plot_map;
const RGBColor WHITE = {255, 255, 255};
} // namespace model
//==============================================================================
// RUN_PLOT controls the logic for making one or many plots
@ -42,7 +47,7 @@ int openmc_plot_geometry()
{
int err;
for (auto pl : plots) {
for (auto pl : model::plots) {
std::stringstream ss;
ss << "Processing plot " << pl.id_ << ": "
<< pl.path_plot_ << "...";
@ -63,11 +68,10 @@ int openmc_plot_geometry()
void
read_plots(pugi::xml_node* plots_node)
{
n_plots = 0;
for (auto node : plots_node->children("plot")) {
Plot pl(node);
plots.push_back(pl);
plot_map[pl.id_] = n_plots++;
model::plots.push_back(pl);
model::plot_map[pl.id_] = model::plots.size() - 1;
}
}
@ -122,7 +126,7 @@ void create_ppm(Plot pl)
p.initialize();
std::copy(xyz, xyz+3, p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = openmc_root_universe;
p.coord[0].universe = model::root_universe;
#pragma omp for
for (int y = 0; y < height; y++) {
@ -155,7 +159,7 @@ Plot::set_id(pugi::xml_node plot_node)
}
// Check to make sure 'id' hasn't been used
if (plot_map.find(id_) != plot_map.end()) {
if (model::plot_map.find(id_) != model::plot_map.end()) {
std::stringstream err_msg;
err_msg << "Two or more plots use the same unique ID: " << id_;
fatal_error(err_msg.str());
@ -243,7 +247,7 @@ Plot::set_bg_color(pugi::xml_node plot_node)
if (check_for_node(plot_node, "background")) {
std::vector<int> bg_rgb = get_node_array<int>(plot_node, "background");
if (PlotType::voxel == type_) {
if (openmc_master) {
if (mpi::master) {
std::stringstream err_msg;
err_msg << "Background color ignored in voxel plot "
<< id_;
@ -361,30 +365,27 @@ Plot::set_default_colors(pugi::xml_node plot_node)
}
if ("cell" == pl_color_by) {
color_by_ = PlotColorBy::cells;
colors_.resize(n_cells);
for (int i = 0; i < n_cells; i++) {
colors_[i] = random_color();
}
colors_.resize(model::cells.size());
} else if("material" == pl_color_by) {
color_by_ = PlotColorBy::mats;
colors_.resize(n_materials);
for (int i = 0; i < materials.size(); i++) {
colors_[i] = random_color();
}
colors_.resize(model::materials.size());
} else {
std::stringstream err_msg;
err_msg << "Unsupported plot color type '" << pl_color_by
<< "' in plot " << id_;
fatal_error(err_msg);
}
for (auto& c : colors_) {
c = random_color();
}
}
void
Plot::set_user_colors(pugi::xml_node plot_node)
{
if (!plot_node.select_nodes("color").empty() && PlotType::voxel == type_) {
if (openmc_master) {
if (mpi::master) {
std::stringstream err_msg;
err_msg << "Color specifications ignored in voxel plot "
<< id_;
@ -412,8 +413,8 @@ Plot::set_user_colors(pugi::xml_node plot_node)
}
// Add RGB
if (PlotColorBy::cells == color_by_) {
if (cell_map.find(col_id) != cell_map.end()) {
col_id = cell_map[col_id];
if (model::cell_map.find(col_id) != model::cell_map.end()) {
col_id = model::cell_map[col_id];
colors_[col_id] = user_rgb;
} else {
std::stringstream err_msg;
@ -422,8 +423,8 @@ Plot::set_user_colors(pugi::xml_node plot_node)
fatal_error(err_msg);
}
} else if (PlotColorBy::mats == color_by_) {
if (material_map.find(col_id) != material_map.end()) {
col_id = material_map[col_id];
if (model::material_map.find(col_id) != model::material_map.end()) {
col_id = model::material_map[col_id];
colors_[col_id] = user_rgb;
} else {
std::stringstream err_msg;
@ -552,7 +553,7 @@ Plot::set_mask(pugi::xml_node plot_node)
if (!mask_nodes.empty()) {
if (PlotType::voxel == type_) {
if (openmc_master) {
if (mpi::master) {
std::stringstream wrn_msg;
wrn_msg << "Mask ignored in voxel plot " << id_;
warning(wrn_msg);
@ -575,8 +576,8 @@ Plot::set_mask(pugi::xml_node plot_node)
// in the cell and material arrays
for (auto& col_id : iarray) {
if (PlotColorBy::cells == color_by_) {
if (cell_map.find(col_id) != cell_map.end()) {
col_id = cell_map[col_id];
if (model::cell_map.find(col_id) != model::cell_map.end()) {
col_id = model::cell_map[col_id];
}
else {
std::stringstream err_msg;
@ -585,8 +586,8 @@ Plot::set_mask(pugi::xml_node plot_node)
fatal_error(err_msg);
}
} else if (PlotColorBy::mats == color_by_) {
if (material_map.find(col_id) != material_map.end()) {
col_id = material_map[col_id];
if (model::material_map.find(col_id) != model::material_map.end()) {
col_id = model::material_map[col_id];
}
else {
std::stringstream err_msg;
@ -660,7 +661,7 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
} else {
if (PlotColorBy::mats == pl.color_by_) {
// Assign color based on material
Cell* c = cells[p.coord[j].cell];
Cell* c = model::cells[p.coord[j].cell];
if (c->type_ == FILL_UNIVERSE) {
// If we stopped on a middle universe level, treat as if not found
rgb = pl.not_found_;
@ -671,12 +672,12 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
id = -1;
} else {
rgb = pl.colors_[p.material - 1];
id = materials[p.material - 1]->id_;
id = model::materials[p.material - 1]->id_;
}
} else if (PlotColorBy::cells == pl.color_by_) {
// Assign color based on cell
rgb = pl.colors_[p.coord[j].cell];
id = cells[p.coord[j].cell]->id_;
id = model::cells[p.coord[j].cell]->id_;
}
} // endif found_cell
}
@ -759,7 +760,7 @@ void draw_mesh_lines(Plot pl, ImageData& data)
width[1] = xyz_ur_plot[1] - xyz_ll_plot[1];
width[2] = xyz_ur_plot[2] - xyz_ll_plot[2];
auto& m = meshes[pl.index_meshlines_mesh_];
auto& m = model::meshes[pl.index_meshlines_mesh_];
int ijk_ll[3], ijk_ur[3];
bool in_mesh;
@ -851,7 +852,7 @@ void create_voxel(Plot pl)
p.initialize();
std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = openmc_root_universe;
p.coord[0].universe = model::root_universe;
// Open binary plot file for writing
std::ofstream of;

4
src/relaxng/materials.rnc
View file

@ -37,7 +37,5 @@ element materials {
}*
}+ &
element cross_sections { xsd:string { maxLength = "255" } }? &
element multipole_library { xsd:string { maxLength = "255" } }?
element cross_sections { xsd:string { maxLength = "255" } }?
}

7
src/relaxng/materials.rng
View file

@ -161,12 +161,5 @@
</data>
</element>
</optional>
<optional>
<element name="multipole_library">
<data type="string">
<param name="maxLength">255</param>
</data>
</element>
</optional>
</interleave>
</element>

19
src/settings.F90
View file

@ -101,7 +101,6 @@ module settings
character(MAX_FILE_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml
character(MAX_FILE_LEN) :: path_multipole ! Path to wmp library
character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point
character(MAX_FILE_LEN) :: path_source_point ! Path to binary source point
character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart
@ -115,7 +114,6 @@ module settings
integer(C_INT), bind(C) :: res_scat_method ! resonance scattering method
real(C_DOUBLE), bind(C) :: res_scat_energy_min
real(C_DOUBLE), bind(C) :: res_scat_energy_max
character(10), allocatable :: res_scat_nuclides(:)
! Is CMFD active
logical(C_BOOL), bind(C) :: cmfd_run
@ -123,21 +121,4 @@ module settings
! No reduction at end of batch
logical(C_BOOL), bind(C) :: reduce_tallies
contains
!===============================================================================
! FREE_MEMORY_SETTINGS deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_settings()
interface
subroutine free_memory_settings_c() bind(C)
end subroutine
end interface
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)
call free_memory_settings_c()
end subroutine free_memory_settings
end module settings

91
src/settings.cpp
View file

@ -17,6 +17,7 @@
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/output.h"
#include "openmc/random_lcg.h"
#include "openmc/simulation.h"
@ -63,7 +64,6 @@ bool dagmc {false};
std::string path_cross_sections;
std::string path_input;
std::string path_multipole;
std::string path_output;
std::string path_particle_restart;
std::string path_source;
@ -73,7 +73,7 @@ std::string path_statepoint;
int32_t index_entropy_mesh {-1};
int32_t index_ufs_mesh {-1};
int32_t index_cmfd_mesh {-1};
int32_t n_batches;
int32_t n_inactive {0};
int32_t gen_per_batch {1};
@ -88,6 +88,7 @@ int n_max_batches;
int res_scat_method {RES_SCAT_ARES};
double res_scat_energy_min {0.01};
double res_scat_energy_max {1000.0};
std::vector<std::string> res_scat_nuclides;
int run_mode {-1};
std::unordered_set<int> sourcepoint_batch;
std::unordered_set<int> statepoint_batch;
@ -150,7 +151,7 @@ void get_run_parameters(pugi::xml_node node_base)
// Preallocate space for keff and entropy by generation
int m = settings::n_max_batches * settings::gen_per_batch;
simulation::k_generation.reserve(m);
entropy.reserve(m);
simulation::entropy.reserve(m);
// Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) {
@ -187,7 +188,7 @@ void read_settings_xml()
using namespace pugi;
// Check if settings.xml exists
std::string filename = std::string(path_input) + "settings.xml";
std::string filename = path_input + "settings.xml";
if (!file_exists(filename)) {
if (run_mode != RUN_MODE_PLOTTING) {
std::stringstream msg;
@ -231,7 +232,7 @@ void read_settings_xml()
// To this point, we haven't displayed any output since we didn't know what
// the verbosity is. Now that we checked for it, show the title if necessary
if (openmc_master) {
if (mpi::master) {
if (verbosity >= 2) title();
}
write_message("Reading settings XML file...", 5);
@ -256,21 +257,6 @@ void read_settings_xml()
path_cross_sections = get_node_value(root, "cross_sections");
}
// Look for deprecated windowed_multipole file in settings.xml
if (run_mode != RUN_MODE_PLOTTING) {
if (check_for_node(root, "multipole_library")) {
warning("Setting multipole_library in settings.xml has been "
"deprecated. The multipole_library is now set in materials.xml and"
" the multipole_library input to materials.xml and the "
"OPENMC_MULTIPOLE_LIBRARY environment variable will take "
"precendent over setting multipole_library in settings.xml.");
path_multipole = get_node_value(root, "multipole_library");
}
if (!ends_with(path_multipole, "/")) {
path_multipole += "/";
}
}
if (!run_CE) {
// Scattering Treatments
if (check_for_node(root, "max_order")) {
@ -413,7 +399,7 @@ void read_settings_xml()
omp_set_num_threads(simulation::n_threads);
}
#else
if (openmc_master) warning("OpenMC was not compiled with OpenMP support; "
if (mpi::master) warning("OpenMC was not compiled with OpenMP support; "
"ignoring number of threads.");
#endif
}
@ -430,17 +416,17 @@ void read_settings_xml()
// Get point to list of <source> elements and make sure there is at least one
for (pugi::xml_node node : root.children("source")) {
external_sources.emplace_back(node);
model::external_sources.emplace_back(node);
}
// If no source specified, default to isotropic point source at origin with Watt spectrum
if (external_sources.empty()) {
if (model::external_sources.empty()) {
SourceDistribution source {
UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})},
UPtrAngle{new Isotropic()},
UPtrDist{new Watt(0.988, 2.249e-6)}
};
external_sources.push_back(std::move(source));
model::external_sources.push_back(std::move(source));
}
// Check if we want to write out source
@ -521,12 +507,12 @@ void read_settings_xml()
// Shannon Entropy mesh
if (check_for_node(root, "entropy_mesh")) {
int temp = std::stoi(get_node_value(root, "entropy_mesh"));
if (mesh_map.find(temp) == mesh_map.end()) {
if (model::mesh_map.find(temp) == model::mesh_map.end()) {
std::stringstream msg;
msg << "Mesh " << temp << " specified for Shannon entropy does not exist.";
fatal_error(msg);
}
index_entropy_mesh = mesh_map.at(temp);
index_entropy_mesh = model::mesh_map.at(temp);
} else if (check_for_node(root, "entropy")) {
warning("Specifying a Shannon entropy mesh via the <entropy> element "
@ -535,18 +521,18 @@ void read_settings_xml()
// Read entropy mesh from <entropy>
auto node_entropy = root.child("entropy");
meshes.emplace_back(new RegularMesh{node_entropy});
model::meshes.emplace_back(new RegularMesh{node_entropy});
// Set entropy mesh index
index_entropy_mesh = meshes.size() - 1;
index_entropy_mesh = model::meshes.size() - 1;
// Assign ID and set mapping
meshes.back()->id_ = 10000;
mesh_map[10000] = index_entropy_mesh;
model::meshes.back()->id_ = 10000;
model::mesh_map[10000] = index_entropy_mesh;
}
if (index_entropy_mesh >= 0) {
auto& m = *meshes[index_entropy_mesh];
auto& m = *model::meshes[index_entropy_mesh];
if (m.shape_.dimension() == 0) {
// If the user did not specify how many mesh cells are to be used in
// each direction, we automatically determine an appropriate number of
@ -566,13 +552,13 @@ void read_settings_xml()
// Uniform fission source weighting mesh
if (check_for_node(root, "ufs_mesh")) {
auto temp = std::stoi(get_node_value(root, "ufs_mesh"));
if (mesh_map.find(temp) == mesh_map.end()) {
if (model::mesh_map.find(temp) == model::mesh_map.end()) {
std::stringstream msg;
msg << "Mesh " << temp << " specified for uniform fission site method "
"does not exist.";
fatal_error(msg);
}
index_ufs_mesh = mesh_map.at(temp);
index_ufs_mesh = model::mesh_map.at(temp);
} else if (check_for_node(root, "uniform_fs")) {
warning("Specifying a UFS mesh via the <uniform_fs> element "
@ -581,14 +567,14 @@ void read_settings_xml()
// Read entropy mesh from <entropy>
auto node_ufs = root.child("uniform_fs");
meshes.emplace_back(new RegularMesh{node_ufs});
model::meshes.emplace_back(new RegularMesh{node_ufs});
// Set entropy mesh index
index_ufs_mesh = meshes.size() - 1;
index_ufs_mesh = model::meshes.size() - 1;
// Assign ID and set mapping
meshes.back()->id_ = 10001;
mesh_map[10001] = index_entropy_mesh;
model::meshes.back()->id_ = 10001;
model::mesh_map[10001] = index_entropy_mesh;
}
if (index_ufs_mesh >= 0) {
@ -749,7 +735,10 @@ void read_settings_xml()
"lower resonance scattering energy bound.");
}
// TODO: Get resonance scattering nuclides
// Get resonance scattering nuclides
if (check_for_node(node_res_scat, "nuclides")) {
res_scat_nuclides = get_node_array<std::string>(node_res_scat, "nuclides");
}
}
// TODO: Get volume calculations
@ -817,22 +806,38 @@ void read_settings_xml()
//==============================================================================
extern "C" {
const char* openmc_path_input() {
bool res_scat_nuclides_empty() {
return settings::res_scat_nuclides.empty();
}
int res_scat_nuclides_size() {
return settings::res_scat_nuclides.size();
}
bool res_scat_nuclides_cmp(int i, const char* name) {
return settings::res_scat_nuclides[i - 1] == name;
}
const char* path_cross_sections_c() {
return settings::path_cross_sections.c_str();
}
const char* path_input_c() {
return settings::path_input.c_str();
}
const char* openmc_path_statepoint() {
const char* path_statepoint_c() {
return settings::path_statepoint.c_str();
}
const char* openmc_path_sourcepoint() {
const char* path_sourcepoint_c() {
return settings::path_sourcepoint.c_str();
}
const char* openmc_path_particle_restart() {
const char* path_particle_restart_c() {
return settings::path_particle_restart.c_str();
}
void free_memory_settings_c() {
void free_memory_settings() {
settings::statepoint_batch.clear();
settings::sourcepoint_batch.clear();
settings::res_scat_nuclides.clear();
}
}

43
src/simulation.cpp
View file

@ -25,10 +25,10 @@ extern "C" bool cmfd_on;
extern "C" void accumulate_tallies();
extern "C" void allocate_banks();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void cmfd_init_batch();
extern "C" void cmfd_tally_init();
extern "C" void configure_tallies();
extern "C" void execute_cmfd();
extern "C" void init_tally_routines();
extern "C" void join_bank_from_threads();
@ -81,7 +81,7 @@ int openmc_simulation_init()
// Allocate array for matching filter bins
#pragma omp parallel
{
filter_matches.resize(n_filters);
simulation::filter_matches.resize(model::tally_filters.size());
}
// Set up tally procedure pointers
@ -92,7 +92,7 @@ int openmc_simulation_init()
allocate_banks();
// Allocate tally results arrays if they're not allocated yet
configure_tallies();
allocate_tally_results();
// Activate the CMFD tallies
cmfd_tally_init();
@ -104,8 +104,9 @@ int openmc_simulation_init()
// will potentially populate k_generation and entropy)
simulation::current_batch = 0;
simulation::k_generation.clear();
entropy.clear();
simulation::entropy.clear();
simulation::need_depletion_rx = false;
openmc_reset();
// If this is a restart run, load the state point data and binary source
// file
@ -139,12 +140,12 @@ int openmc_simulation_finalize()
if (!simulation::initialized) return 0;
// Stop active batch timer and start finalization timer
time_active.stop();
time_finalize.start();
simulation::time_active.stop();
simulation::time_finalize.start();
#pragma omp parallel
{
filter_matches.clear();
simulation::filter_matches.clear();
}
// Deallocate Fortran variables, set tallies to inactive
@ -166,8 +167,8 @@ int openmc_simulation_finalize()
}
// Stop timers and show timing statistics
time_finalize.stop();
time_total.stop();
simulation::time_finalize.stop();
simulation::time_total.stop();
if (mpi::master) {
if (settings::verbosity >= 6) print_runtime();
if (settings::verbosity >= 4) print_results();
@ -200,7 +201,7 @@ int openmc_next_batch(int* status)
initialize_generation();
// Start timer for transport
time_transport.start();
simulation::time_transport.start();
// ====================================================================
// LOOP OVER PARTICLES
@ -218,7 +219,7 @@ int openmc_next_batch(int* status)
}
// Accumulate time for transport
time_transport.stop();
simulation::time_transport.stop();
finalize_generation();
}
@ -305,10 +306,10 @@ void initialize_batch()
// Manage active/inactive timers and activate tallies if necessary.
if (first_inactive) {
time_inactive.start();
simulation::time_inactive.start();
} else if (first_active) {
time_inactive.stop();
time_active.start();
simulation::time_inactive.stop();
simulation::time_active.start();
for (int i = 1; i <= n_tallies; ++i) {
// TODO: change one-based index
openmc_tally_set_active(i, true);
@ -327,9 +328,9 @@ void initialize_batch()
void finalize_batch()
{
// Reduce tallies onto master process and accumulate
time_tallies.start();
simulation::time_tallies.start();
accumulate_tallies();
time_tallies.stop();
simulation::time_tallies.stop();
// Reset global tally results
if (simulation::current_batch <= settings::n_inactive) {
@ -384,13 +385,13 @@ void initialize_generation()
{
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Reset number of fission bank sites
n_bank = 0;
simulation::n_bank = 0;
// Count source sites if using uniform fission source weighting
if (settings::ufs_on) ufs_count_sites();
// Store current value of tracklength k
keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
simulation::keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
}
}
@ -564,4 +565,10 @@ extern "C" void k_generation_clear() { simulation::k_generation.clear(); }
extern "C" void k_generation_reserve(int i) { simulation::k_generation.reserve(i); }
extern "C" int64_t work_index(int rank) { return simulation::work_index[rank]; }
// This function was moved here to get around a bug on macOS whereby an invalid
// pointer is returned for the threadprivate filter_matches
extern "C" FilterMatch* filter_match_pointer(int indx) {
return &simulation::filter_matches[indx];
}
} // namespace openmc

34
src/source.cpp
View file

@ -27,8 +27,12 @@ namespace openmc {
// Global variables
//==============================================================================
namespace model {
std::vector<SourceDistribution> external_sources;
}
//==============================================================================
// SourceDistribution implementation
//==============================================================================
@ -167,9 +171,9 @@ Bank SourceDistribution::sample() const
if (space_box) {
if (space_box->only_fissionable()) {
// Determine material
auto c = cells[cell_index - 1];
auto c = model::cells[cell_index - 1];
int32_t mat_index = c->material_[instance];
auto m = materials[mat_index];
auto m = model::materials[mat_index];
if (mat_index == MATERIAL_VOID) {
found = false;
@ -207,10 +211,10 @@ Bank SourceDistribution::sample() const
auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
if (energy_ptr) {
auto energies = xt::adapt(energy_ptr->x());
if (xt::any(energies > energy_max[p-1])) {
if (xt::any(energies > data::energy_max[p-1])) {
fatal_error("Source energy above range of energies of at least "
"one cross section table");
} else if (xt::any(energies < energy_min[p-1])) {
} else if (xt::any(energies < data::energy_min[p-1])) {
fatal_error("Source energy below range of energies of at least "
"one cross section table");
}
@ -221,7 +225,7 @@ Bank SourceDistribution::sample() const
site.E = energy_->sample();
// Resample if energy falls outside minimum or maximum particle energy
if (site.E < energy_max[p-1] && site.E > energy_min[p-1]) break;
if (site.E < data::energy_max[p-1] && site.E > data::energy_min[p-1]) break;
}
// Set delayed group
@ -299,28 +303,28 @@ Bank sample_external_source()
// Determine total source strength
double total_strength = 0.0;
for (auto& s : external_sources)
for (auto& s : model::external_sources)
total_strength += s.strength();
// Sample from among multiple source distributions
int i = 0;
if (external_sources.size() > 1) {
if (model::external_sources.size() > 1) {
double xi = prn()*total_strength;
double c = 0.0;
for (; i < external_sources.size(); ++i) {
c += external_sources[i].strength();
for (; i < model::external_sources.size(); ++i) {
c += model::external_sources[i].strength();
if (xi < c) break;
}
}
// Sample source site from i-th source distribution
Bank site {external_sources[i].sample()};
Bank site {model::external_sources[i].sample()};
// If running in MG, convert site % E to group
if (!settings::run_CE) {
site.E = lower_bound_index(rev_energy_bins.begin(), rev_energy_bins.end(),
site.E);
site.E = num_energy_groups - site.E;
site.E = lower_bound_index(data::rev_energy_bins.begin(),
data::rev_energy_bins.end(), site.E);
site.E = data::num_energy_groups - site.E;
}
// Set the random number generator back to the tracking stream.
@ -335,13 +339,13 @@ Bank sample_external_source()
extern "C" void free_memory_source()
{
external_sources.clear();
model::external_sources.clear();
}
extern "C" double total_source_strength()
{
double strength = 0.0;
for (const auto& s : external_sources) {
for (const auto& s : model::external_sources) {
strength += s.strength();
}
return strength;

4
src/state_point.F90
View file

@ -55,7 +55,7 @@ contains
!===============================================================================
function openmc_statepoint_write(filename, write_source) result(err) bind(C)
type(C_PTR), intent(in), optional :: filename
type(C_PTR), value :: filename
logical(C_BOOL), intent(in), optional :: write_source
integer(C_INT) :: err
@ -91,7 +91,7 @@ contains
err = 0
! Set the filename
if (present(filename)) then
if (c_associated(filename)) then
call c_f_pointer(filename, string, [MAX_FILE_LEN])
filename_ = to_f_string(string)
else

8
src/state_point.cpp
View file

@ -113,7 +113,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
#else
if (openmc_master) {
if (mpi::master) {
// Create dataset big enough to hold all source sites
hsize_t dims[] {static_cast<hsize_t>(settings::n_particles)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
@ -320,11 +320,11 @@ void restart_set_keff()
{
if (simulation::restart_batch > settings::n_inactive) {
for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
k_sum[0] += simulation::k_generation[i];
k_sum[1] += std::pow(simulation::k_generation[i], 2);
simulation::k_sum[0] += simulation::k_generation[i];
simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch*n_realizations;
simulation::keff = k_sum[0] / n;
simulation::keff = simulation::k_sum[0] / n;
} else {
simulation::keff = simulation::k_generation.back();
}

40
src/string_functions.cpp
View file

@ -1,40 +0,0 @@
#include "openmc/string_functions.h"
#include <sstream>
namespace openmc {
std::string& strtrim(std::string& s)
{
const char* t = " \t\n\r\f\v";
s.erase(s.find_last_not_of(t) + 1);
s.erase(0, s.find_first_not_of(t));
return s;
}
char* strtrim(char* c_str)
{
std::string std_str;
std_str.assign(c_str);
strtrim(std_str);
int length = std_str.copy(c_str, std_str.size());
c_str[length] = '\0';
return c_str;
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
int word_count(std::string const& str)
{
std::stringstream stream(str);
std::string dum;
int count = 0;
while (stream >> dum) {count++;}
return count;
}
} // namespace openmc

78
src/string_utils.cpp Normal file
View file

@ -0,0 +1,78 @@
#include "openmc/string_utils.h"
#include <algorithm> // for equal
#include <cctype> // for tolower, isspace
#include <sstream>
namespace openmc {
std::string& strtrim(std::string& s)
{
const char* t = " \t\n\r\f\v";
s.erase(s.find_last_not_of(t) + 1);
s.erase(0, s.find_first_not_of(t));
return s;
}
char* strtrim(char* c_str)
{
std::string std_str;
std_str.assign(c_str);
strtrim(std_str);
int length = std_str.copy(c_str, std_str.size());
c_str[length] = '\0';
return c_str;
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
}
int word_count(std::string const& str)
{
std::stringstream stream(str);
std::string dum;
int count = 0;
while (stream >> dum) {count++;}
return count;
}
std::vector<std::string> split(const std::string& in)
{
std::vector<std::string> out;
for (int i = 0; i < in.size(); ) {
// Increment i until we find a non-whitespace character.
if (std::isspace(in[i])) {
i++;
} else {
// Find the next whitespace character at j.
int j = i + 1;
while (j < in.size() && std::isspace(in[j]) == 0) {j++;}
// Push-back everything between i and j.
out.push_back(in.substr(i, j-i));
i = j + 1; // j is whitespace so leapfrog to j+1
}
}
return out;
}
bool ends_with(const std::string& value, const std::string& ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
bool starts_with(const std::string& value, const std::string& beginning)
{
if (beginning.size() > value.size()) return false;
return std::equal(beginning.begin(), beginning.end(), value.begin());
}
} // namespace openmc

18
src/summary.cpp
View file

@ -10,7 +10,7 @@ extern "C" void
write_geometry(hid_t file_id) {
auto geom_group = create_group(file_id, "geometry");
#ifdef DAGMC
if (settings::dagmc) {
write_attribute(geom_group, "dagmc", 1);
@ -18,25 +18,25 @@ write_geometry(hid_t file_id) {
}
#endif
write_attribute(geom_group, "n_cells", cells.size());
write_attribute(geom_group, "n_surfaces", surfaces.size());
write_attribute(geom_group, "n_universes", universes.size());
write_attribute(geom_group, "n_lattices", lattices.size());
write_attribute(geom_group, "n_cells", model::cells.size());
write_attribute(geom_group, "n_surfaces", model::surfaces.size());
write_attribute(geom_group, "n_universes", model::universes.size());
write_attribute(geom_group, "n_lattices", model::lattices.size());
auto cells_group = create_group(geom_group, "cells");
for (Cell* c : cells) c->to_hdf5(cells_group);
for (Cell* c : model::cells) c->to_hdf5(cells_group);
close_group(cells_group);
auto surfaces_group = create_group(geom_group, "surfaces");
for (Surface* surf : surfaces) surf->to_hdf5(surfaces_group);
for (Surface* surf : model::surfaces) surf->to_hdf5(surfaces_group);
close_group(surfaces_group);
auto universes_group = create_group(geom_group, "universes");
for (Universe* u : universes) u->to_hdf5(universes_group);
for (Universe* u : model::universes) u->to_hdf5(universes_group);
close_group(universes_group);
auto lattices_group = create_group(geom_group, "lattices");
for (Lattice* lat : lattices) lat->to_hdf5(lattices_group);
for (Lattice* lat : model::lattices) lat->to_hdf5(lattices_group);
close_group(lattices_group);
close_group(geom_group);

75
src/surface.cpp
View file

@ -8,7 +8,7 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xml_interface.h"
#include "openmc/string_functions.h"
#include "openmc/string_utils.h"
namespace openmc {
@ -25,12 +25,13 @@ extern "C" const int BC_PERIODIC {3};
// Global variables
//==============================================================================
int32_t n_surfaces;
namespace model {
std::vector<Surface*> surfaces;
std::map<int, int> surface_map;
} // namespace model
//==============================================================================
// Helper functions for reading the "coeffs" node of an XML surface element
//==============================================================================
@ -1063,13 +1064,14 @@ extern "C" void
read_surfaces(pugi::xml_node* node)
{
// Count the number of surfaces.
for (pugi::xml_node surf_node: node->children("surface")) {n_surfaces++;}
int n_surfaces = 0;
for (pugi::xml_node surf_node : node->children("surface")) {n_surfaces++;}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
// Loop over XML surface elements and populate the array.
surfaces.reserve(n_surfaces);
model::surfaces.reserve(n_surfaces);
{
pugi::xml_node surf_node;
int i_surf;
@ -1078,40 +1080,40 @@ read_surfaces(pugi::xml_node* node)
std::string surf_type = get_node_value(surf_node, "type", true, true);
if (surf_type == "x-plane") {
surfaces.push_back(new SurfaceXPlane(surf_node));
model::surfaces.push_back(new SurfaceXPlane(surf_node));
} else if (surf_type == "y-plane") {
surfaces.push_back(new SurfaceYPlane(surf_node));
model::surfaces.push_back(new SurfaceYPlane(surf_node));
} else if (surf_type == "z-plane") {
surfaces.push_back(new SurfaceZPlane(surf_node));
model::surfaces.push_back(new SurfaceZPlane(surf_node));
} else if (surf_type == "plane") {
surfaces.push_back(new SurfacePlane(surf_node));
model::surfaces.push_back(new SurfacePlane(surf_node));
} else if (surf_type == "x-cylinder") {
surfaces.push_back(new SurfaceXCylinder(surf_node));
model::surfaces.push_back(new SurfaceXCylinder(surf_node));
} else if (surf_type == "y-cylinder") {
surfaces.push_back(new SurfaceYCylinder(surf_node));
model::surfaces.push_back(new SurfaceYCylinder(surf_node));
} else if (surf_type == "z-cylinder") {
surfaces.push_back(new SurfaceZCylinder(surf_node));
model::surfaces.push_back(new SurfaceZCylinder(surf_node));
} else if (surf_type == "sphere") {
surfaces.push_back(new SurfaceSphere(surf_node));
model::surfaces.push_back(new SurfaceSphere(surf_node));
} else if (surf_type == "x-cone") {
surfaces.push_back(new SurfaceXCone(surf_node));
model::surfaces.push_back(new SurfaceXCone(surf_node));
} else if (surf_type == "y-cone") {
surfaces.push_back(new SurfaceYCone(surf_node));
model::surfaces.push_back(new SurfaceYCone(surf_node));
} else if (surf_type == "z-cone") {
surfaces.push_back(new SurfaceZCone(surf_node));
model::surfaces.push_back(new SurfaceZCone(surf_node));
} else if (surf_type == "quadric") {
surfaces.push_back(new SurfaceQuadric(surf_node));
model::surfaces.push_back(new SurfaceQuadric(surf_node));
} else {
std::stringstream err_msg;
@ -1122,11 +1124,11 @@ read_surfaces(pugi::xml_node* node)
}
// Fill the surface map.
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
int id = surfaces[i_surf]->id_;
auto in_map = surface_map.find(id);
if (in_map == surface_map.end()) {
surface_map[id] = i_surf;
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
int id = model::surfaces[i_surf]->id_;
auto in_map = model::surface_map.find(id);
if (in_map == model::surface_map.end()) {
model::surface_map[id] = i_surf;
} else {
std::stringstream err_msg;
err_msg << "Two or more surfaces use the same unique ID: " << id;
@ -1138,10 +1140,10 @@ read_surfaces(pugi::xml_node* node)
double xmin {INFTY}, xmax {-INFTY}, ymin {INFTY}, ymax {-INFTY},
zmin {INFTY}, zmax {-INFTY};
int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax;
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
if (surfaces[i_surf]->bc_ == BC_PERIODIC) {
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = surfaces[i_surf];
Surface* surf_base = model::surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Make sure this surface inherits from PeriodicSurface.
@ -1183,10 +1185,10 @@ read_surfaces(pugi::xml_node* node)
}
// Set i_periodic for periodic BC surfaces.
for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
if (surfaces[i_surf]->bc_ == BC_PERIODIC) {
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = surfaces[i_surf];
Surface* surf_base = model::surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Also try downcasting to the SurfacePlane type (which must be handled
@ -1216,7 +1218,7 @@ read_surfaces(pugi::xml_node* node)
}
} else {
// Convert the surface id to an index.
surf->i_periodic_ = surface_map[surf->i_periodic_];
surf->i_periodic_ = model::surface_map[surf->i_periodic_];
}
} else {
// This is a SurfacePlane. We won't try to find it's partner if the
@ -1228,12 +1230,12 @@ read_surfaces(pugi::xml_node* node)
fatal_error(err_msg);
} else {
// Convert the surface id to an index.
surf->i_periodic_ = surface_map[surf->i_periodic_];
surf->i_periodic_ = model::surface_map[surf->i_periodic_];
}
}
// Make sure the opposite surface is also periodic.
if (surfaces[surf->i_periodic_]->bc_ != BC_PERIODIC) {
if (model::surfaces[surf->i_periodic_]->bc_ != BC_PERIODIC) {
std::stringstream err_msg;
err_msg << "Could not find matching surface for periodic boundary "
"condition on surface " << surf->id_;
@ -1248,7 +1250,7 @@ read_surfaces(pugi::xml_node* node)
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind) {return surfaces[surf_ind];}
Surface* surface_pointer(int surf_ind) {return model::surfaces[surf_ind];}
int surface_id(Surface* surf) {return surf->id_;}
@ -1288,11 +1290,12 @@ extern "C" {
void free_memory_surfaces_c()
{
for (Surface* surf : surfaces) {delete surf;}
surfaces.clear();
n_surfaces = 0;
surface_map.clear();
for (Surface* surf : model::surfaces) {delete surf;}
model::surfaces.clear();
model::surface_map.clear();
}
int surfaces_size() { return model::surfaces.size(); }
}
} // namespace openmc

9
src/surface_header.F90
View file

@ -78,10 +78,15 @@ module surface_header
end type Surface
integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces
type(Surface), allocatable, target :: surfaces(:)
interface
function surfaces_size() result(sz) bind(C)
import C_INT
integer(C_INT) :: sz
end function
end interface
contains
pure function surface_id(this) result(id)

Some files were not shown because too many files have changed in this diff Show more