crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' of https://github.com/mit-crpg/openmc into new-update

Browse source 
merging commit 1004
This commit is contained in:
Jose Salcedo Perez 2018-05-12 01:10:26 +00:00
commit 2452747a9b
158 changed files with 8888 additions and 10229 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.gitignore vendored
View file

@ -99,4 +99,4 @@ examples/jupyter/plots
.tox/
.python-version
.coverage
htmlcov
htmlcov

36
CMakeLists.txt
View file

@ -5,14 +5,10 @@ project(openmc Fortran C CXX)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR}/include)
# Set module path
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
# Make sure Fortran module directory is included when building
include_directories(${CMAKE_BINARY_DIR}/include)
#===============================================================================
# Architecture specific definitions
#===============================================================================
@ -73,7 +69,7 @@ if(NOT DEFINED HDF5_PREFER_PARALLEL)
endif()
endif()
find_package(HDF5 COMPONENTS Fortran_HL)
find_package(HDF5 COMPONENTS HL)
if(NOT HDF5_FOUND)
message(FATAL_ERROR "Could not find HDF5")
endif()
@ -436,18 +432,22 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/trigger_header.F90
)
set(LIBOPENMC_CXX_SRC
src/error.h
src/hdf5_interface.h
src/initialize.cpp
src/finalize.cpp
src/hdf5_interface.cpp
src/message_passing.cpp
src/plot.cpp
src/random_lcg.cpp
src/random_lcg.h
src/simulation.cpp
src/state_point.cpp
src/surface.cpp
src/surface.h
src/xml_interface.h
src/pugixml/pugixml.cpp
src/pugixml/pugixml.hpp)
src/xml_interface.cpp
src/pugixml/pugixml.cpp)
add_library(libopenmc SHARED ${LIBOPENMC_FORTRAN_SRC} ${LIBOPENMC_CXX_SRC})
set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc)
add_executable(${program} src/main.F90)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
PUBLIC_HEADER include/openmc.h)
add_executable(${program} src/main.cpp)
#===============================================================================
# Add compiler/linker flags
@ -456,11 +456,11 @@ add_executable(${program} src/main.F90)
set_property(TARGET ${program} libopenmc pugixml_fortran
PROPERTY LINKER_LANGUAGE Fortran)
target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS})
target_include_directories(libopenmc PUBLIC include ${HDF5_INCLUDE_DIRS})
# The executable and the faddeeva package use only one language. They can be
# set via target_compile_options which accepts a list.
target_compile_options(${program} PUBLIC ${f90flags})
target_compile_options(${program} PUBLIC ${cxxflags})
target_compile_options(faddeeva PRIVATE ${cflags})
# The libopenmc library has both F90 and C++ so the compile flags must be set
@ -501,7 +501,9 @@ add_custom_command(TARGET libopenmc POST_BUILD
install(TARGETS ${program} libopenmc
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib)
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include
)
install(DIRECTORY src/relaxng DESTINATION share/openmc)
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)

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

@ -46,10 +46,13 @@ Type Definitions
Functions
---------
.. c:function:: void openmc_calculate_volumes()
.. c:function:: int openmc_calculate_volumes()
Run a stochastic volume calculation
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
Get the fill for a cell
@ -192,11 +195,14 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_finalize()
.. c:function:: int openmc_finalize()
Finalize a simulation
.. c:function:: void openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance)
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance)
Determine the ID of the cell/material containing a given point
@ -207,6 +213,8 @@ Functions
occurs, the ID is -1.
:param int32_t* instance: If a cell is repeated in the geometry, the instance
of the cell that was found and zero otherwise.
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_cell_index(int32_t id, int32_t* index)
@ -266,17 +274,24 @@ Functions
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_hard_reset()
.. c:function:: int openmc_hard_reset()
Reset tallies, timers, and pseudo-random number generator state
.. c:function:: void openmc_init(const int* intracomm)
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_init(int argc, char** argv, const void* intracomm)
Initialize OpenMC
:param int argc: Number of command-line arguments (including command)
:param char** argv: Command-line arguments
:param intracomm: MPI intracommunicator. If MPI is not being used, a null
pointer should be passed.
:type intracomm: const int*
:type intracomm: const void*
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_load_nuclide(char name[])
@ -394,26 +409,41 @@ Functions
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_plot_geometry()
.. c:function:: int openmc_plot_geometry()
Run plotting mode.
.. c:function:: void openmc_reset()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_reset()
Resets all tally scores
.. c:function:: void openmc_run()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_run()
Run a simulation
.. c:function:: void openmc_simulation_finalize()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_simulation_finalize()
Finalize a simulation.
.. c:function:: void openmc_simulation_init()
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_simulation_init()
Initialize a simulation. Must be called after openmc_init().
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_source_bank(struct Bank** ptr, int64_t* n)
Return a pointer to the source bank array.
@ -433,13 +463,15 @@ Functions
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_statepoint_write(const char filename[])
.. c:function:: int openmc_statepoint_write(const char filename[])
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[]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_tally_get_id(int32_t index, int32_t* id)

87
docs/source/io_formats/settings.rst
View file

@ -93,29 +93,13 @@ or ``multi-group``.
*Default*: continuous-energy
---------------------
``<entropy>`` Element
---------------------
--------------------------
``<entropy_mesh>`` Element
--------------------------
The ``<entropy>`` element describes a mesh that is used for calculating Shannon
entropy. This mesh should cover all possible fissionable materials in the
problem. It has the following attributes/sub-elements:
:dimension:
The number of mesh cells in the x, y, and z directions, respectively.
*Default*: If this tag is not present, the number of mesh cells is
automatically determined by the code.
:lower_left:
The Cartesian coordinates of the lower-left corner of the mesh.
*Default*: None
:upper_right:
The Cartesian coordinates of the upper-right corner of the mesh.
*Default*: None
The ``<entropy_mesh>`` element indicates the ID of a mesh that is to be used for
calculating Shannon entropy. The mesh should cover all possible fissionable
materials in the problem and is specified using a :ref:`mesh_element`.
-----------------------------------
``<generations_per_batch>`` Element
@ -199,6 +183,36 @@ then, OpenMC will only use up to the :math:`P_1` data.
.. note:: This element is not used in the continuous-energy
:ref:`energy_mode`.
.. _mesh_element:
------------------
``<mesh>`` Element
------------------
The ``<mesh>`` element describes a mesh that is used either for calculating
Shannon entropy, applying the uniform fission site method, or in tallies. For
Shannon entropy meshes, the mesh should cover all possible fissionable materials
in the problem. It has the following attributes/sub-elements:
:id:
A unique integer that is used to identify the mesh.
:dimension:
The number of mesh cells in the x, y, and z directions, respectively.
*Default*: If this tag is not present, the number of mesh cells is
automatically determined by the code.
:lower_left:
The Cartesian coordinates of the lower-left corner of the mesh.
*Default*: None
:upper_right:
The Cartesian coordinates of the upper-right corner of the mesh.
*Default*: None
-----------------------
``<no_reduce>`` Element
-----------------------
@ -765,30 +779,15 @@ has the following attributes/sub-elements:
------------------------
``<uniform_fs>`` Element
``<ufs_mesh>`` Element
------------------------
The ``<uniform_fs>`` element describes a mesh that is used for re-weighting
source sites at every generation based on the uniform fission site methodology
described in Kelly et al., "MC21 Analysis of the Nuclear Energy Agency Monte
Carlo Performance Benchmark Problem," Proceedings of *Physor 2012*, Knoxville,
TN (2012). This mesh should cover all possible fissionable materials in the
problem. It has the following attributes/sub-elements:
:dimension:
The number of mesh cells in the x, y, and z directions, respectively.
*Default*: None
:lower_left:
The Cartesian coordinates of the lower-left corner of the mesh.
*Default*: None
:upper_right:
The Cartesian coordinates of the upper-right corner of the mesh.
*Default*: None
The ``<ufs_mesh>`` element indicates the ID of a mesh that is used for
re-weighting source sites at every generation based on the uniform fission site
methodology described in Kelly et al., "MC21 Analysis of the Nuclear Energy
Agency Monte Carlo Performance Benchmark Problem," Proceedings of *Physor 2012*,
Knoxville, TN (2012). The mesh should cover all possible fissionable materials
in the problem and is specified using a :ref:`mesh_element`.
.. _verbosity:

16
docs/source/io_formats/source.rst
View file

@ -8,13 +8,13 @@ Normally, source data is stored in a state point file. However, it is possible
to request that the source be written separately, in which case the format used
is that documented here.
**/filetype** (*char[]*)
**/**
String indicating the type of file.
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
**/source_bank** (Compound type)
Source bank information for each particle. The compound type has fields
``wgt``, ``xyz``, ``uvw``, ``E``, and ``delayed_group``, which
represent the weight, position, direction, energy, energy group, and
delayed_group of the source particle, respectively.
:Datasets:
- **source_bank** (Compound type) -- Source bank information for each
particle. The compound type has fields ``wgt``, ``xyz``, ``uvw``,
``E``, and ``delayed_group``, which represent the weight, position,
direction, energy, energy group, and delayed_group of the source
particle, respectively.

9
docs/source/io_formats/statepoint.rst
View file

@ -133,15 +133,8 @@ The current version of the statepoint file format is 17.0.
- **derivative** (*int*) -- ID of the derivative applied to the
tally.
- **n_score_bins** (*int*) -- Number of scoring bins for a single
nuclide. In general, this can be greater than the number of
user-specified scores since each score might have multiple scoring
bins, e.g., scatter-PN.
nuclide.
- **score_bins** (*char[][]*) -- Values of specified scores.
- **n_user_scores** (*int*) -- Number of scores without accounting
for those added by expansions, e.g. scatter-PN.
- **moment_orders** (*char[][]*) -- Tallying moment orders for
Legendre and spherical harmonic tally expansions (e.g., 'P2',
'Y1,2', etc.).
- **results** (*double[][][2]*) -- Accumulated sum and sum-of-squares
for each bin of the i-th tally. The first dimension represents
combinations of filter bins, the second dimensions represents

16
docs/source/io_formats/summary.rst
View file

@ -4,7 +4,7 @@
Summary File Format
===================
The current version of the summary file format is 5.0.
The current version of the summary file format is 6.0.
**/**
@ -104,8 +104,13 @@ The current version of the summary file format is 5.0.
- **atom_density** (*double[]*) -- Total atom density of the material
in atom/b-cm.
- **nuclides** (*char[][]*) -- Array of nuclides present in the
material, e.g., 'U235'.
material, e.g., 'U235'. This data set is only present if nuclides
are used.
- **nuclide_densities** (*double[]*) -- Atom density of each nuclide.
This data set is only present if 'nuclides' data set is present.
- **macroscopics** (*char[][]*) -- Array of macroscopic data sets
present in the material. This dataset is only present if
macroscopic data sets are used in multi-group mode.
- **sab_names** (*char[][]*) -- Names of
S(:math:`\alpha,\beta`) tables assigned to the material.
@ -116,6 +121,13 @@ The current version of the summary file format is 5.0.
:Datasets: - **names** (*char[][]*) -- Names of nuclides.
- **awrs** (*float[]*) -- Atomic weight ratio of each nuclide.
**/macroscopics/**
:Attributes: - **n_macroscopics** (*int*) -- Number of macroscopic data sets
in the problem.
:Datasets: - **names** (*char[][]*) -- Names of the macroscopic data sets.
**/tallies/tally <uid>/**
:Datasets: - **name** (*char[]*) -- Name of the tally.

3
docs/source/pythonapi/capi.rst
View file

@ -44,5 +44,8 @@ Classes
EnergyFilter
MaterialFilter
Material
Mesh
MeshFilter
MeshSurfaceFilter
Nuclide
Tally

8
docs/source/usersguide/install.rst
View file

@ -141,16 +141,10 @@ Prerequisites
recommend that your HDF5 installation be built with parallel I/O
features. An example of configuring HDF5_ is listed below::
FC=mpifort ./configure --enable-fortran --enable-parallel
FC=mpifort ./configure --enable-parallel
You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial.
.. important::
If you are building HDF5 version 1.8.x or earlier, you must include
``--enable-fortran2003`` as well when configuring HDF5 or else OpenMC
will not be able to compile.
.. admonition:: Optional
:class: note

63
docs/source/usersguide/tallies.rst
View file

@ -21,7 +21,12 @@ region of phase space, as in:
Thus, to specify a tally, we need to specify what regions of phase space should
be included when deciding whether to score an event as well as what the scoring
function (:math:`f` in the above equation) should be used. The regions of phase
space are called *filters* and the scoring functions are simply called *scores*.
space are generally called *filters* and the scoring functions are simply
called *scores*.
The only cases when filters do not correspond directly with the regions of
phase space are when expansion functions are applied in the integrand, such as
for Legendre expansions of the scattering kernel.
-------
Filters
@ -69,10 +74,9 @@ Scores
------
To specify the scoring functions, a list of strings needs to be given to the
:attr:`Tally.scores` attribute. You can score the flux ('flux'), a reaction rate
('total', 'fission', etc.), or even scattering moments (e.g., 'scatter-P3'). For
example, to tally the elastic scattering rate and the fission neutron
production, you'd assign::
:attr:`Tally.scores` attribute. You can score the flux ('flux'), or a reaction
rate ('total', 'fission', etc.). For example, to tally the elastic scattering
rate and the fission neutron production, you'd assign::
tally.scores = ['elastic', 'nu-fission']
@ -98,12 +102,6 @@ The following tables show all valid scores:
+======================+===================================================+
|flux |Total flux. |
+----------------------+---------------------------------------------------+
|flux-YN |Spherical harmonic expansion of the direction of |
| |motion :math:`\left(\Omega\right)` of the total |
| |flux. This score will tally all of the harmonic |
| |moments of order 0 to N. N must be between 0 and |
| |10. |
+----------------------+---------------------------------------------------+
.. table:: **Reaction scores: units are reactions per source particle.**
@ -118,43 +116,10 @@ The following tables show all valid scores:
+----------------------+---------------------------------------------------+
|fission |Total fission reaction rate. |
+----------------------+---------------------------------------------------+
|scatter |Total scattering rate. Can also be identified with |
| |the "scatter-0" response type. |
+----------------------+---------------------------------------------------+
|scatter-N |Tally the N\ :sup:`th` \ scattering moment, where N|
| |is the Legendre expansion order of the change in |
| |particle angle :math:`\left(\mu\right)`. N must be |
| |between 0 and 10. As an example, tallying the 2\ |
| |:sup:`nd` \ scattering moment would be specified as|
| |``<scores>scatter-2</scores>``. |
+----------------------+---------------------------------------------------+
|scatter-PN |Tally all of the scattering moments from order 0 to|
| |N, where N is the Legendre expansion order of the |
| |change in particle angle |
| |:math:`\left(\mu\right)`. That is, "scatter-P1" is |
| |equivalent to requesting tallies of "scatter-0" and|
| |"scatter-1". Like for "scatter-N", N must be |
| |between 0 and 10. As an example, tallying up to the|
| |2\ :sup:`nd` \ scattering moment would be specified|
| |as ``<scores> scatter-P2 </scores>``. |
+----------------------+---------------------------------------------------+
|scatter-YN |"scatter-YN" is similar to "scatter-PN" except an |
| |additional expansion is performed for the incoming |
| |particle direction :math:`\left(\Omega\right)` |
| |using the real spherical harmonics. This is useful|
| |for performing angular flux moment weighting of the|
| |scattering moments. Like "scatter-PN", "scatter-YN"|
| |will tally all of the moments from order 0 to N; N |
| |again must be between 0 and 10. |
|scatter |Total scattering rate. |
+----------------------+---------------------------------------------------+
|total |Total reaction rate. |
+----------------------+---------------------------------------------------+
|total-YN |The total reaction rate expanded via spherical |
| |harmonics about the direction of motion of the |
| |neutron, :math:`\Omega`. This score will tally all |
| |of the harmonic moments of order 0 to N. N must be|
| |between 0 and 10. |
+----------------------+---------------------------------------------------+
|(n,2nd) |(n,2nd) reaction rate. |
+----------------------+---------------------------------------------------+
|(n,2n) |(n,2n) reaction rate. |
@ -248,10 +213,10 @@ The following tables show all valid scores:
+----------------------+---------------------------------------------------+
|nu-fission |Total production of neutrons due to fission. |
+----------------------+---------------------------------------------------+
|nu-scatter, |These scores are similar in functionality to their |
|nu-scatter-N, |``scatter*`` equivalents except the total |
|nu-scatter-PN, |production of neutrons due to scattering is scored |
|nu-scatter-YN |vice simply the scattering rate. This accounts for |
|nu-scatter, |This score is similar in functionality to the |
| |``scatter`` score except the total production of |
| |neutrons due to scattering is scored vice simply |
| |the scattering rate. This accounts for |
| |multiplicity from (n,2n), (n,3n), and (n,4n) |
| |reactions. |
+----------------------+---------------------------------------------------+

411
examples/jupyter/mg-mode-part-ii.ipynb
View file

File diff suppressed because one or more lines are too long

309
examples/jupyter/mg-mode-part-iii.ipynb
View file

File diff suppressed because one or more lines are too long

919
examples/jupyter/mgxs-part-iii.ipynb
View file

File diff suppressed because one or more lines are too long

90
include/openmc.h
View file

@ -16,7 +16,7 @@ extern "C" {
int delayed_group;
};
void openmc_calculate_voumes();
int openmc_calculate_volumes();
int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
int openmc_cell_get_id(int32_t index, int32_t* id);
int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
@ -27,21 +27,27 @@ extern "C" {
int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_sources(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_filter_get_id(int32_t index, int32_t* id);
int openmc_filter_get_type(int32_t index, const char** type);
int openmc_filter_set_id(int32_t index, int32_t id);
void openmc_finalize();
int openmc_filter_set_type(int32_t index, const char* type);
int openmc_finalize();
int openmc_find(double* xyz, int rtype, int32_t* id, int32_t* instance);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
void openmc_get_filter_next_id(int32_t* id);
int openmc_get_keff(double k_combined[]);
int openmc_get_material_index(int32_t id, int32_t* index);
int openmc_get_mesh_index(int32_t id, int32_t* index);
int openmc_get_nuclide_index(const char name[], int* index);
int64_t openmc_get_seed();
int openmc_get_tally_index(int32_t id, int32_t* index);
void openmc_hard_reset();
void openmc_init(const int* intracomm);
int openmc_hard_reset();
int openmc_init(int argc, char* argv[], const void* intracomm);
int openmc_init_f(const int* intracomm);
int openmc_legendre_filter_get_order(int32_t index, int* order);
int openmc_legendre_filter_set_order(int32_t index, int order);
int openmc_load_nuclide(char name[]);
@ -53,15 +59,25 @@ extern "C" {
int openmc_material_set_id(int32_t index, int32_t id);
int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n);
int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins);
int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_mesh_get_id(int32_t index, int32_t* id);
int openmc_mesh_get_dimension(int32_t index, int** id, int* n);
int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
int openmc_mesh_set_id(int32_t index, int32_t id);
int openmc_mesh_set_dimension(int32_t index, int n, const int* dims);
int openmc_mesh_set_params(int32_t index, const double* ll, const double* ur, const double* width, int n);
int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_next_batch();
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, char** name);
void openmc_plot_geometry();
void openmc_reset();
void openmc_run();
void openmc_simulation_finalize();
void openmc_simulation_init();
int openmc_particle_restart();
int openmc_plot_geometry();
int openmc_reset();
int openmc_run();
void openmc_set_seed(int64_t new_seed);
int openmc_simulation_finalize();
int openmc_simulation_init();
int openmc_source_bank(struct Bank** ptr, int64_t* n);
int openmc_source_set_strength(int32_t index, double strength);
int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
@ -73,17 +89,20 @@ extern "C" {
int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
void openmc_statepoint_write(const char filename[]);
int openmc_statepoint_write(const char filename[]);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n);
int openmc_tally_get_n_realizations(int32_t index, int32_t* n);
int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
int openmc_tally_get_scores(int32_t index, int** scores, int* n);
int openmc_tally_results(int32_t index, double** ptr, int shape_[3]);
int openmc_tally_set_active(int32_t index, bool active);
int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices);
int openmc_tally_set_id(int32_t index, int32_t id);
int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
int openmc_tally_set_scores(int32_t index, int n, const int* scores);
int openmc_tally_set_scores(int32_t index, int n, const char** scores);
int openmc_tally_set_type(int32_t index, const char* type);
int openmc_zernike_filter_get_order(int32_t index, int* order);
int openmc_zernike_filter_get_params(int32_t index, double* x, double* y, double* r);
int openmc_zernike_filter_set_order(int32_t index, int order);
@ -91,22 +110,22 @@ extern "C" {
const double* y, const double* r);
// Error codes
extern int E_UNASSIGNED;
extern int E_ALLOCATE;
extern int E_OUT_OF_BOUNDS;
extern int E_INVALID_SIZE;
extern int E_INVALID_ARGUMENT;
extern int E_INVALID_TYPE;
extern int E_INVALID_ID;
extern int E_GEOMETRY;
extern int E_DATA;
extern int E_PHYSICS;
extern int E_WARNING;
extern int OPENMC_E_UNASSIGNED;
extern int OPENMC_E_ALLOCATE;
extern int OPENMC_E_OUT_OF_BOUNDS;
extern int OPENMC_E_INVALID_SIZE;
extern int OPENMC_E_INVALID_ARGUMENT;
extern int OPENMC_E_INVALID_TYPE;
extern int OPENMC_E_INVALID_ID;
extern int OPENMC_E_GEOMETRY;
extern int OPENMC_E_DATA;
extern int OPENMC_E_PHYSICS;
extern int OPENMC_E_WARNING;
// Global variables
extern char openmc_err_msg[256];
extern double keff;
extern double keff_std;
extern double openmc_keff;
extern double openmc_keff_std;
extern int32_t n_batches;
extern int32_t n_cells;
extern int32_t n_filters;
@ -122,9 +141,24 @@ extern "C" {
extern int32_t n_surfaces;
extern int32_t n_tallies;
extern int32_t n_universes;
extern int run_mode;
extern bool simulation_initialized;
extern int verbosity;
extern int openmc_run_mode;
extern bool openmc_simulation_initialized;
extern int openmc_verbosity;
// Variables that are shared by necessity (can be removed from public header
// later)
extern bool openmc_master;
extern int openmc_n_procs;
extern int openmc_n_threads;
extern int openmc_rank;
extern int64_t openmc_work;
// Run modes
constexpr int RUN_MODE_FIXEDSOURCE {1};
constexpr int RUN_MODE_EIGENVALUE {2};
constexpr int RUN_MODE_PLOTTING {3};
constexpr int RUN_MODE_PARTICLE {4};
constexpr int RUN_MODE_VOLUME {5};
#ifdef __cplusplus
}

1
openmc/capi/__init__.py
View file

@ -43,6 +43,7 @@ from .core import *
from .nuclide import *
from .material import *
from .cell import *
from .mesh import *
from .filter import *
from .tally import *
from .settings import settings

18
openmc/capi/cell.py
View file

@ -5,9 +5,10 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError, InvalidIDError
from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .error import _error_handler
from .material import Material
__all__ = ['Cell', 'cells']
@ -44,7 +45,7 @@ class Cell(_FortranObjectWithID):
This class exposes a cell that is stored internally in the OpenMC
library. To obtain a view of a cell with a given ID, use the
:data:`openmc.capi.nuclides` mapping.
:data:`openmc.capi.cells` mapping.
Parameters
----------
@ -115,14 +116,15 @@ class Cell(_FortranObjectWithID):
def fill(self, fill):
if isinstance(fill, Iterable):
n = len(fill)
indices = (c_int*n)(*(m._index for m in fill))
_dll.openmc_cell_set_fill(self._index, 1, 1, indices)
indices = (c_int32*n)(*(m._index if m is not None else -1
for m in fill))
_dll.openmc_cell_set_fill(self._index, 1, n, indices)
elif isinstance(fill, Material):
materials = [fill]
indices = (c_int*1)(fill._index)
indices = (c_int32*1)(fill._index)
_dll.openmc_cell_set_fill(self._index, 1, 1, indices)
elif fill is None:
indices = (c_int32*1)(-1)
_dll.openmc_cell_set_fill(self._index, 1, 1, indices)
else:
raise NotImplementedError
def set_temperature(self, T, instance=None):
"""Set the temperature of a cell

105
openmc/capi/core.py
View file

@ -1,13 +1,14 @@
from contextlib import contextmanager
from ctypes import (CDLL, c_int, c_int32, c_int64, c_double, c_char_p,
POINTER, Structure)
from ctypes import (CDLL, c_int, c_int32, c_int64, c_double, c_char_p, c_char,
POINTER, Structure, c_void_p, create_string_buffer)
from warnings import warn
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError
from . import _dll
from .error import _error_handler, AllocationError
from .error import _error_handler
import openmc.capi
@ -19,29 +20,41 @@ class _Bank(Structure):
('delayed_group', c_int)]
_dll.openmc_calculate_volumes.restype = None
_dll.openmc_finalize.restype = None
_dll.openmc_calculate_volumes.restype = c_int
_dll.openmc_calculate_volumes.errcheck = _error_handler
_dll.openmc_finalize.restype = c_int
_dll.openmc_finalize.errcheck = _error_handler
_dll.openmc_find.argtypes = [POINTER(c_double*3), c_int, POINTER(c_int32),
POINTER(c_int32)]
_dll.openmc_find.restype = c_int
_dll.openmc_find.errcheck = _error_handler
_dll.openmc_hard_reset.restype = None
_dll.openmc_init.argtypes = [POINTER(c_int)]
_dll.openmc_init.restype = None
_dll.openmc_hard_reset.restype = c_int
_dll.openmc_hard_reset.errcheck = _error_handler
_dll.openmc_init.argtypes = [c_int, POINTER(POINTER(c_char)), c_void_p]
_dll.openmc_init.restype = c_int
_dll.openmc_init.errcheck = _error_handler
_dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
_dll.openmc_get_keff.restype = c_int
_dll.openmc_get_keff.errcheck = _error_handler
_dll.openmc_next_batch.argtypes = [POINTER(c_int)]
_dll.openmc_next_batch.restype = c_int
_dll.openmc_plot_geometry.restype = None
_dll.openmc_run.restype = None
_dll.openmc_reset.restype = None
_dll.openmc_next_batch.errcheck = _error_handler
_dll.openmc_plot_geometry.restype = c_int
_dll.openmc_plot_geometry.restype = _error_handler
_dll.openmc_run.restype = c_int
_dll.openmc_run.errcheck = _error_handler
_dll.openmc_reset.restype = c_int
_dll.openmc_reset.errcheck = _error_handler
_dll.openmc_source_bank.argtypes = [POINTER(POINTER(_Bank)), POINTER(c_int64)]
_dll.openmc_source_bank.restype = c_int
_dll.openmc_source_bank.errcheck = _error_handler
_dll.openmc_simulation_init.restype = None
_dll.openmc_simulation_finalize.restype = None
_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)]
_dll.openmc_statepoint_write.restype = None
_dll.openmc_statepoint_write.restype = c_int
_dll.openmc_statepoint_write.errcheck = _error_handler
def calculate_volumes():
@ -102,25 +115,42 @@ def hard_reset():
_dll.openmc_hard_reset()
def init(intracomm=None):
def init(args=None, intracomm=None):
"""Initialize OpenMC
Parameters
----------
args : list of str
Command-line arguments
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
"""
if intracomm is not None:
# If an mpi4py communicator was passed, convert it to an integer to
# be passed to openmc_init
try:
intracomm = intracomm.py2f()
except AttributeError:
pass
_dll.openmc_init(c_int(intracomm))
if args is not None:
args = ['openmc'] + list(args)
argc = len(args)
# Create the argv array. Note that it is actually expected to be of
# length argc + 1 with the final item being a null pointer.
argv = (POINTER(c_char) * (argc + 1))()
for i, arg in enumerate(args):
argv[i] = create_string_buffer(arg.encode())
else:
_dll.openmc_init(None)
argc = 0
argv = None
if intracomm is not None:
# If an mpi4py communicator was passed, convert it to void* to be passed
# to openmc_init
try:
from mpi4py import MPI
except ImportError:
intracomm = None
else:
address = MPI._addressof(intracomm)
intracomm = c_void_p(address)
_dll.openmc_init(argc, argv, intracomm)
def iter_batches():
@ -147,13 +177,13 @@ def iter_batches():
"""
while True:
# Run next batch
retval = next_batch()
status = next_batch()
# Provide opportunity for user to perform action between batches
yield
# End the iteration
if retval < 0:
if status != 0:
break
@ -174,18 +204,25 @@ def keff():
return tuple(k)
else:
# Otherwise, return the tracklength estimator
mean = c_double.in_dll(_dll, 'keff').value
std_dev = c_double.in_dll(_dll, 'keff_std').value if n > 1 else np.inf
mean = c_double.in_dll(_dll, 'openmc_keff').value
std_dev = c_double.in_dll(_dll, 'openmc_keff_std').value \
if n > 1 else np.inf
return (mean, std_dev)
def next_batch():
"""Run next batch."""
retval = _dll.openmc_next_batch()
if retval == -3:
raise AllocationError('Simulation has not been initialized. You must call '
'openmc.capi.simulation_init() first.')
return retval
"""Run next batch.
Returns
-------
int
Status after running a batch (0=normal, 1=reached maximum number of
batches, 2=tally triggers reached)
"""
status = c_int()
_dll.openmc_next_batch(status)
return status.value
def plot_geometry():

73
openmc/capi/error.py
View file

@ -1,45 +1,10 @@
from ctypes import c_int, c_char
from warnings import warn
import openmc.exceptions as exc
from . import _dll
class OpenMCError(Exception):
"""Root exception class for OpenMC."""
class GeometryError(OpenMCError):
"""Geometry-related error"""
class InvalidIDError(OpenMCError):
"""Use of an ID that is invalid."""
class AllocationError(OpenMCError):
"""Error related to memory allocation."""
class OutOfBoundsError(OpenMCError):
"""Index in array out of bounds."""
class DataError(OpenMCError):
"""Error relating to nuclear data."""
class PhysicsError(OpenMCError):
"""Error relating to performing physics."""
class InvalidArgumentError(OpenMCError):
"""Argument passed was invalid."""
class InvalidTypeError(OpenMCError):
"""Tried to perform an operation on the wrong type."""
def _error_handler(err, func, args):
"""Raise exception according to error code."""
@ -52,23 +17,23 @@ def _error_handler(err, func, args):
msg = errmsg.value.decode()
# Raise exception type corresponding to error code
if err == errcode('e_allocate'):
raise AllocationError(msg)
elif err == errcode('e_out_of_bounds'):
raise OutOfBoundsError(msg)
elif err == errcode('e_invalid_argument'):
raise InvalidArgumentError(msg)
elif err == errcode('e_invalid_type'):
raise InvalidTypeError(msg)
if err == errcode('e_invalid_id'):
raise InvalidIDError(msg)
elif err == errcode('e_geometry'):
raise GeometryError(msg)
elif err == errcode('e_data'):
raise DataError(msg)
elif err == errcode('e_physics'):
raise PhysicsError(msg)
elif err == errcode('e_warning'):
if err == errcode('OPENMC_E_ALLOCATE'):
raise exc.AllocationError(msg)
elif err == errcode('OPENMC_E_OUT_OF_BOUNDS'):
raise exc.OutOfBoundsError(msg)
elif err == errcode('OPENMC_E_INVALID_ARGUMENT'):
raise exc.InvalidArgumentError(msg)
elif err == errcode('OPENMC_E_INVALID_TYPE'):
raise exc.InvalidTypeError(msg)
if err == errcode('OPENMC_E_INVALID_ID'):
raise exc.InvalidIDError(msg)
elif err == errcode('OPENMC_E_GEOMETRY'):
raise exc.GeometryError(msg)
elif err == errcode('OPENMC_E_DATA'):
raise exc.DataError(msg)
elif err == errcode('OPENMC_E_PHYSICS'):
raise exc.PhysicsError(msg)
elif err == errcode('OPENMC_E_WARNING'):
warn(msg)
elif err < 0:
raise OpenMCError("Unknown error encountered (code {}).".format(err))
raise exc.OpenMCError("Unknown error encountered (code {}).".format(err))

42
openmc/capi/filter.py
View file

@ -6,17 +6,19 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError, InvalidIDError
from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .error import _error_handler
from .material import Material
from .mesh import Mesh
__all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
'CellbornFilter', 'CellfromFilter', 'DistribcellFilter',
'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter',
'EnergyFunctionFilter', 'MaterialFilter', 'MeshFilter',
'MuFilter', 'PolarFilter', 'SurfaceFilter',
'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SurfaceFilter',
'UniverseFilter', 'filters']
# Tally functions
@ -52,9 +54,15 @@ _dll.openmc_material_filter_get_bins.errcheck = _error_handler
_dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_int32)]
_dll.openmc_material_filter_set_bins.restype = c_int
_dll.openmc_material_filter_set_bins.errcheck = _error_handler
_dll.openmc_mesh_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_mesh_filter_get_mesh.restype = c_int
_dll.openmc_mesh_filter_get_mesh.errcheck = _error_handler
_dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
_dll.openmc_mesh_filter_set_mesh.restype = c_int
_dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
@ -190,10 +198,40 @@ class MaterialFilter(Filter):
class MeshFilter(Filter):
filter_type = 'mesh'
def __init__(self, mesh=None, uid=None, new=True, index=None):
super().__init__(uid, new, index)
if mesh is not None:
self.mesh = mesh
@property
def mesh(self):
index_mesh = c_int32()
_dll.openmc_mesh_filter_get_mesh(self._index, index_mesh)
return Mesh(index=index_mesh.value)
@mesh.setter
def mesh(self, mesh):
_dll.openmc_mesh_filter_set_mesh(self._index, mesh._index)
class MeshSurfaceFilter(Filter):
filter_type = 'meshsurface'
def __init__(self, mesh=None, uid=None, new=True, index=None):
super().__init__(uid, new, index)
if mesh is not None:
self.mesh = mesh
@property
def mesh(self):
index_mesh = c_int32()
_dll.openmc_meshsurface_filter_get_mesh(self._index, index_mesh)
return Mesh(index=index_mesh.value)
@mesh.setter
def mesh(self, mesh):
_dll.openmc_meshsurface_filter_set_mesh(self._index, mesh._index)
class MuFilter(Filter):
filter_type = 'mu'

6
openmc/capi/material.py
View file

@ -5,9 +5,10 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError, InvalidIDError
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .error import _error_handler
__all__ = ['Material', 'materials']
@ -89,6 +90,9 @@ class Material(_FortranObjectWithID):
index = index.value
else:
index = mapping[uid]._index
elif index == -1:
# Special value indicates void material
return None
if index not in cls.__instances:
instance = super(Material, cls).__new__(cls)

183
openmc/capi/mesh.py Normal file
View file

@ -0,0 +1,183 @@
from collections.abc import Mapping, Iterable
from ctypes import c_int, c_int32, c_double, POINTER
from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError, InvalidIDError
from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
from .material import Material
__all__ = ['Mesh', 'meshes']
# Mesh functions
_dll.openmc_extend_meshes.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
_dll.openmc_extend_meshes.restype = c_int
_dll.openmc_extend_meshes.errcheck = _error_handler
_dll.openmc_mesh_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_mesh_get_id.restype = c_int
_dll.openmc_mesh_get_id.errcheck = _error_handler
_dll.openmc_mesh_get_dimension.argtypes = [c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_mesh_get_dimension.restype = c_int
_dll.openmc_mesh_get_dimension.errcheck = _error_handler
_dll.openmc_mesh_get_params.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(POINTER(c_double)),
POINTER(POINTER(c_double)), POINTER(c_int)]
_dll.openmc_mesh_get_params.restype = c_int
_dll.openmc_mesh_get_params.errcheck = _error_handler
_dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_mesh_set_id.restype = c_int
_dll.openmc_mesh_set_id.errcheck = _error_handler
_dll.openmc_mesh_set_dimension.argtypes = [c_int32, c_int, POINTER(c_int)]
_dll.openmc_mesh_set_dimension.restype = c_int
_dll.openmc_mesh_set_dimension.errcheck = _error_handler
_dll.openmc_mesh_set_params.argtypes = [
c_int32, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double)]
_dll.openmc_mesh_set_params.restype = c_int
_dll.openmc_mesh_set_params.errcheck = _error_handler
_dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_mesh_index.restype = c_int
_dll.openmc_get_mesh_index.errcheck = _error_handler
class Mesh(_FortranObjectWithID):
"""Mesh stored internally.
This class exposes a mesh that is stored internally in the OpenMC
library. To obtain a view of a mesh with a given ID, use the
:data:`openmc.capi.meshes` mapping.
Parameters
----------
index : int
Index in the `meshes` array.
Attributes
----------
id : int
ID of the mesh
dimension : iterable of int
The number of mesh cells in each direction.
lower_left : numpy.ndarray
The lower-left corner of the structured mesh. If only two coordinate are
given, it is assumed that the mesh is an x-y mesh.
upper_right : numpy.ndarray
The upper-right corner of the structrued mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
width : numpy.ndarray
The width of mesh cells in each direction.
"""
__instances = WeakValueDictionary()
def __new__(cls, uid=None, new=True, index=None):
mapping = meshes
if index is None:
if new:
# Determine ID to assign
if uid is None:
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A mesh with ID={} has already '
'been allocated.'.format(uid))
index = c_int32()
_dll.openmc_extend_meshes(1, index, None)
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
mesh_id = c_int32()
_dll.openmc_mesh_get_id(self._index, mesh_id)
return mesh_id.value
@id.setter
def id(self, mesh_id):
_dll.openmc_mesh_set_id(self._index, mesh_id)
@property
def dimension(self):
dims = POINTER(c_int)()
n = c_int()
_dll.openmc_mesh_get_dimension(self._index, dims, n)
return tuple(as_array(dims, (n.value,)))
@dimension.setter
def dimension(self, dimension):
n = len(dimension)
dimension = (c_int*n)(*dimension)
_dll.openmc_mesh_set_dimension(self._index, n, dimension)
@property
def lower_left(self):
return self._get_parameters()[0]
@property
def upper_right(self):
return self._get_parameters()[1]
@property
def width(self):
return self._get_parameters()[2]
def _get_parameters(self):
ll = POINTER(c_double)()
ur = POINTER(c_double)()
w = POINTER(c_double)()
n = c_int()
_dll.openmc_mesh_get_params(self._index, ll, ur, w, n)
return (
as_array(ll, (n.value,)),
as_array(ur, (n.value,)),
as_array(w, (n.value,))
)
def set_parameters(self, lower_left=None, upper_right=None, width=None):
if lower_left is not None:
n = len(lower_left)
lower_left = (c_double*n)(*lower_left)
if upper_right is not None:
n = len(upper_right)
upper_right = (c_double*n)(*upper_right)
if width is not None:
n = len(width)
width = (c_double*n)(*width)
_dll.openmc_mesh_set_params(self._index, n, lower_left, upper_right, width)
class _MeshMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
try:
_dll.openmc_get_mesh_index(key, index)
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return Mesh(index=index.value)
def __iter__(self):
for i in range(len(self)):
yield Mesh(index=i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_meshes').value
def __repr__(self):
return repr(dict(self))
meshes = _MeshMapping()

3
openmc/capi/nuclide.py
View file

@ -5,9 +5,10 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import DataError, AllocationError
from . import _dll
from .core import _FortranObject
from .error import _error_handler, DataError, AllocationError
from .error import _error_handler
__all__ = ['Nuclide', 'nuclides', 'load_nuclide']

6
openmc/capi/settings.py
View file

@ -20,11 +20,11 @@ class _Settings(object):
generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch')
inactive = _DLLGlobal(c_int32, 'n_inactive')
particles = _DLLGlobal(c_int64, 'n_particles')
verbosity = _DLLGlobal(c_int, 'verbosity')
verbosity = _DLLGlobal(c_int, 'openmc_verbosity')
@property
def run_mode(self):
i = c_int.in_dll(_dll, 'run_mode').value
i = c_int.in_dll(_dll, 'openmc_run_mode').value
try:
return _RUN_MODES[i]
except KeyError:
@ -32,7 +32,7 @@ class _Settings(object):
@run_mode.setter
def run_mode(self, mode):
current_idx = c_int.in_dll(_dll, 'run_mode')
current_idx = c_int.in_dll(_dll, 'openmc_run_mode')
for idx, mode_value in _RUN_MODES.items():
if mode_value == mode:
current_idx.value = idx

29
openmc/capi/tally.py
View file

@ -1,15 +1,16 @@
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER
from ctypes import c_int, c_int32, c_double, c_char_p, c_bool, POINTER
from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
import scipy.stats
from openmc.exceptions import AllocationError, InvalidIDError
from openmc.data.reaction import REACTION_NAME
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .error import _error_handler
from .filter import _get_filter
@ -25,6 +26,9 @@ _dll.openmc_get_tally_index.errcheck = _error_handler
_dll.openmc_global_tallies.argtypes = [POINTER(POINTER(c_double))]
_dll.openmc_global_tallies.restype = c_int
_dll.openmc_global_tallies.errcheck = _error_handler
_dll.openmc_tally_get_active.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_active.restype = c_int
_dll.openmc_tally_get_active.errcheck = _error_handler
_dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_id.restype = c_int
_dll.openmc_tally_get_id.errcheck = _error_handler
@ -47,6 +51,9 @@ _dll.openmc_tally_results.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)]
_dll.openmc_tally_results.restype = c_int
_dll.openmc_tally_results.errcheck = _error_handler
_dll.openmc_tally_set_active.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_active.restype = c_int
_dll.openmc_tally_set_active.errcheck = _error_handler
_dll.openmc_tally_set_filters.argtypes = [c_int32, c_int, POINTER(c_int32)]
_dll.openmc_tally_set_filters.restype = c_int
_dll.openmc_tally_set_filters.errcheck = _error_handler
@ -66,10 +73,10 @@ _dll.openmc_tally_set_type.errcheck = _error_handler
_SCORES = {
-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter',
-9: 'absorption', -10: 'fission', -11: 'nu-fission', -12: 'kappa-fission',
-13: 'current', -18: 'events', -19: 'delayed-nu-fission',
-20: 'prompt-nu-fission', -21: 'inverse-velocity', -22: 'fission-q-prompt',
-23: 'fission-q-recoverable', -24: 'decay-rate'
-5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
-9: 'current', -10: 'events', -11: 'delayed-nu-fission',
-12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-15: 'fission-q-recoverable', -16: 'decay-rate'
}
@ -177,6 +184,16 @@ class Tally(_FortranObjectWithID):
return cls.__instances[index]
@property
def active(self):
active = c_bool()
_dll.openmc_tally_get_active(self._index, active)
return active.value
@active.setter
def active(self, active):
_dll.openmc_tally_set_active(self._index, active)
@property
def id(self):
tally_id = c_int32()

9
openmc/data/library.py
View file

@ -23,9 +23,14 @@ class DataLibrary(EqualityMixin):
def __init__(self):
self.libraries = []
def get_by_material(self, value):
def get_by_material(self, name):
"""Return the library dictionary containing a given material.
Parameters
----------
name : str
Name of material, e.g. 'Am241'
Returns
-------
library : dict or None
@ -34,7 +39,7 @@ class DataLibrary(EqualityMixin):
"""
for library in self.libraries:
if value in library['materials']:
if name in library['materials']:
return library
return None

2
openmc/deplete/operator.py
View file

@ -268,7 +268,7 @@ class Operator(TransportOperator):
# Initialize OpenMC library
comm.barrier()
openmc.capi.init(comm)
openmc.capi.init(intracomm=comm)
# Generate tallies in memory
self._generate_tallies()

38
openmc/exceptions.py Normal file
View file

@ -0,0 +1,38 @@
class OpenMCError(Exception):
"""Root exception class for OpenMC."""
class GeometryError(OpenMCError):
"""Geometry-related error"""
class InvalidIDError(OpenMCError):
"""Use of an ID that is invalid."""
class AllocationError(OpenMCError):
"""Error related to memory allocation."""
class OutOfBoundsError(OpenMCError):
"""Index in array out of bounds."""
class DataError(OpenMCError):
"""Error relating to nuclear data."""
class PhysicsError(OpenMCError):
"""Error relating to performing physics."""
class InvalidArgumentError(OpenMCError):
"""Argument passed was invalid."""
class InvalidTypeError(OpenMCError):
"""Tried to perform an operation on the wrong type."""
class SetupError(OpenMCError):
"""Error while setting up a problem."""

4
openmc/filter.py
View file

@ -104,10 +104,8 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
return False
elif len(self.bins) != len(other.bins):
return False
elif not np.allclose(self.bins, other.bins):
return False
else:
return True
return np.allclose(self.bins, other.bins)
def __ne__(self, other):
return not self == other

20
openmc/filter_expansion.py
View file

@ -10,10 +10,17 @@ from . import Filter
class ExpansionFilter(Filter):
"""Abstract filter class for functional expansions."""
def __init__(self, order, filter_id=None):
self.order = order
self.id = filter_id
def __eq__(self, other):
if type(self) is not type(other):
return False
else:
return self.bins == other.bins
@property
def order(self):
return self._order
@ -318,15 +325,15 @@ class ZernikeFilter(ExpansionFilter):
This filter allows scores to be multiplied by Zernike polynomials of the
particle's position normalized to a given unit circle, up to a
user-specified order. The Zernike polynomials follow the definition by `Noll
<https://doi.org/10.1364/JOSA.66.000207>`_ and are defined as
user-specified order. The standard Zernike polynomials follow the definition by
Born and Wolf, *Principles of Optics* and are defined as
.. math::
Z_n^m(\rho, \theta) = \sqrt{2n + 2} R_n^m(\rho) \cos (m\theta), \quad m > 0
Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0
Z_n^{m}(\rho, \theta) = \sqrt{2n + 2} R_n^{m}(\rho) \sin (m\theta), \quad m < 0
Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0
Z_n^{m}(\rho, \theta) = \sqrt{n + 1} R_n^{m}(\rho), \quad m = 0
Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0
where the radial polynomials are
@ -335,7 +342,8 @@ class ZernikeFilter(ExpansionFilter):
\frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
is exactly :math:`\pi` for each polynomial.
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is
2 if :math:`m` equals 0 and 1 otherwise.
Specifying a filter with order N tallies moments for all :math:`n` from 0 to
N and each value of :math:`m`. The ordering of the Zernike polynomial

20
openmc/material.py
View file

@ -278,8 +278,8 @@ class Material(IDManagerMixin):
name = group['name'].value.decode() if 'name' in group else ''
density = group['atom_density'].value
nuc_densities = group['nuclide_densities'][...]
nuclides = group['nuclides'].value
if 'nuclide_densities' in group:
nuc_densities = group['nuclide_densities'][...]
# Create the Material
material = cls(mat_id, name)
@ -295,10 +295,18 @@ class Material(IDManagerMixin):
# Set the Material's density to atom/b-cm as used by OpenMC
material.set_density(density=density, units='atom/b-cm')
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
name = fullname.decode().strip()
material.add_nuclide(name, percent=density, percent_type='ao')
if 'nuclides' in group:
nuclides = group['nuclides'].value
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
name = fullname.decode().strip()
material.add_nuclide(name, percent=density, percent_type='ao')
if 'macroscopics' in group:
macroscopics = group['macroscopics'].value
# Add all macroscopics to the Material
for fullname in macroscopics:
name = fullname.decode().strip()
material.add_macroscopic(name)
return material

68
openmc/mesh.py
View file

@ -182,41 +182,6 @@ class Mesh(IDManagerMixin):
return mesh
def cell_generator(self):
"""Generator function to traverse through every [i,j,k] index of the
mesh
For example the following code:
.. code-block:: python
for mesh_index in mymesh.cell_generator():
print(mesh_index)
will produce the following output for a 3-D 2x2x2 mesh in mymesh::
[1, 1, 1]
[2, 1, 1]
[1, 2, 1]
[2, 2, 1]
...
"""
if len(self.dimension) == 1:
for x in range(self.dimension[0]):
yield [x + 1, 1, 1]
elif len(self.dimension) == 2:
for y in range(self.dimension[1]):
for x in range(self.dimension[0]):
yield [x + 1, y + 1, 1]
else:
for z in range(self.dimension[2]):
for y in range(self.dimension[1]):
for x in range(self.dimension[0]):
yield [x + 1, y + 1, z + 1]
def to_xml_element(self):
"""Return XML representation of the mesh
@ -280,12 +245,14 @@ class Mesh(IDManagerMixin):
cv.check_value('bc', entry, ['transmission', 'vacuum',
'reflective', 'periodic'])
n_dim = len(self.dimension)
# Build the cell which will contain the lattice
xplanes = [openmc.XPlane(x0=self.lower_left[0],
boundary_type=bc[0]),
openmc.XPlane(x0=self.upper_right[0],
boundary_type=bc[1])]
if len(self.dimension) == 1:
if n_dim == 1:
yplanes = [openmc.YPlane(y0=-1e10, boundary_type='reflective'),
openmc.YPlane(y0=1e10, boundary_type='reflective')]
else:
@ -294,7 +261,7 @@ class Mesh(IDManagerMixin):
openmc.YPlane(y0=self.upper_right[1],
boundary_type=bc[3])]
if len(self.dimension) <= 2:
if n_dim <= 2:
# Would prefer to have the z ranges be the max supported float, but
# these values are apparently different between python and Fortran.
# Choosing a safe and sane default.
@ -314,12 +281,12 @@ class Mesh(IDManagerMixin):
(+yplanes[0] & -yplanes[1]) &
(+zplanes[0] & -zplanes[1]))
# Build the universes which will be used for each of the [i,j,k]
# Build the universes which will be used for each of the (i,j,k)
# locations within the mesh.
# We will concurrently build cells to assign to these universes
cells = []
universes = []
for [i, j, k] in self.cell_generator():
for index in self.indices:
cells.append(openmc.Cell())
universes.append(openmc.Universe())
universes[-1].add_cell(cells[-1])
@ -329,7 +296,24 @@ class Mesh(IDManagerMixin):
# Assign the universe and rotate to match the indexing expected for
# the lattice
lattice.universes = np.rot90(np.reshape(universes, self.dimension))
if n_dim == 1:
universe_array = np.array([universes])
elif n_dim == 2:
universe_array = np.empty(self.dimension, dtype=openmc.Universe)
i = 0
for y in range(self.dimension[1] - 1, -1, -1):
for x in range(self.dimension[0]):
universe_array[y][x] = universes[i]
i += 1
else:
universe_array = np.empty(self.dimension, dtype=openmc.Universe)
i = 0
for z in range(self.dimension[2]):
for y in range(self.dimension[1] - 1, -1, -1):
for x in range(self.dimension[0]):
universe_array[z][y][x] = universes[i]
i += 1
lattice.universes = universe_array
if self.width is not None:
lattice.pitch = self.width
@ -337,9 +321,9 @@ class Mesh(IDManagerMixin):
dx = ((self.upper_right[0] - self.lower_left[0]) /
self.dimension[0])
if len(self.dimension) == 1:
if n_dim == 1:
lattice.pitch = [dx]
elif len(self.dimension) == 2:
elif n_dim == 2:
dy = ((self.upper_right[1] - self.lower_left[1]) /
self.dimension[1])
lattice.pitch = [dx, dy]

88
openmc/mgxs/library.py
View file

@ -1220,9 +1220,8 @@ class Library(object):
xs_type = 'macro'
# Initialize file
mgxs_file = openmc.MGXSLibrary(self.energy_groups,
num_delayed_groups=\
self.num_delayed_groups)
mgxs_file = openmc.MGXSLibrary(
self.energy_groups, num_delayed_groups=self.num_delayed_groups)
if self.domain_type == 'mesh':
# Create the xsdata objects and add to the mgxs_file
@ -1231,7 +1230,7 @@ class Library(object):
if self.by_nuclide:
raise NotImplementedError("Mesh domains do not currently "
"support nuclidic tallies")
for subdomain in domain.cell_generator():
for subdomain in domain.indices:
# Build & add metadata to XSdata object
if xsdata_names is None:
xsdata_name = 'set' + str(i + 1)
@ -1346,7 +1345,7 @@ class Library(object):
geometry.root_universe = root
materials = openmc.Materials()
for i, subdomain in enumerate(self.domains[0].cell_generator()):
for i, subdomain in enumerate(self.domains[0].indices):
xsdata = mgxs_file.xsdatas[i]
# Build the macroscopic and assign it to the cell of
@ -1401,24 +1400,16 @@ class Library(object):
The rules to check include:
- Either total or transport should be present.
- Either total or transport must be present.
- Both can be available if one wants, but we should
use whatever corresponds to Library.correction (if P0: transport)
- Absorption and total (or transport) are required.
- Absorption is required.
- A nu-fission cross section and chi values are not required as a
fixed source problem could be the target.
- Fission and kappa-fission are not required as they are only
needed to support tallies the user may wish to request.
- A nu-scatter matrix is required.
- Having a multiplicity matrix is preferred.
- Having both nu-scatter (of any order) and scatter
(at least isotropic) matrices is the second choice.
- If only nu-scatter, need total (not transport), to
be used in adjusting absorption
(i.e., reduced_abs = tot - nuscatt)
See also
--------
@ -1428,36 +1419,51 @@ class Library(object):
"""
error_flag = False
# if correction is 'P0', then transport must be provided
# otherwise total must be provided
if self.correction == 'P0':
if ('transport' not in self.mgxs_types and
'nu-transport' not in self.mgxs_types):
error_flag = True
warn('If the "correction" parameter is "P0", then a '
'"transport" or "nu-transport" MGXS type is required.')
else:
if 'total' not in self.mgxs_types:
error_flag = True
warn('If the "correction" parameter is None, then a '
'"total" MGXS type is required.')
# Check consistency of "nu-transport" and "nu-scatter"
if 'nu-transport' in self.mgxs_types:
if not ('nu-scatter matrix' in self.mgxs_types or
'consistent nu-scatter matrix' in self.mgxs_types):
error_flag = True
warn('If a "nu-transport" MGXS type is used then a '
'"nu-scatter matrix" or "consistent nu-scatter matrix" '
'must also be used.')
elif 'transport' in self.mgxs_types:
if not ('scatter matrix' in self.mgxs_types or
'consistent scatter matrix' in self.mgxs_types):
error_flag = True
warn('If a "transport" MGXS type is used then a '
'"scatter matrix" or "consistent scatter matrix" '
'must also be used.')
# Make sure there is some kind of a scattering matrix data
if 'nu-scatter matrix' not in self.mgxs_types and \
'consistent nu-scatter matrix' not in self.mgxs_types and \
'scatter matrix' not in self.mgxs_types and \
'consistent scatter matrix' not in self.mgxs_types:
error_flag = True
warn('A "nu-scatter matrix", "consistent nu-scatter matrix", '
'"scatter matrix", or "consistent scatter matrix" MGXS '
'type is required.')
# Ensure absorption is present
if 'absorption' not in self.mgxs_types:
error_flag = True
warn('An "absorption" MGXS type is required but not provided.')
# Ensure nu-scattering matrix is required
if 'nu-scatter matrix' not in self.mgxs_types and \
'consistent nu-scatter matrix' not in self.mgxs_types:
error_flag = True
warn('A "nu-scatter matrix" MGXS type is required but not provided.')
else:
# Ok, now see the status of scatter and/or multiplicity
if 'scatter matrix' not in self.mgxs_types or \
'consistent scatter matrix' not in self.mgxs_types and \
'multiplicity matrix' not in self.mgxs_types:
# We dont have data needed for multiplicity matrix, therefore
# we need total, and not transport.
if 'total' not in self.mgxs_types:
error_flag = True
warn('A "total" MGXS type is required if a '
'scattering matrix is not provided.')
# Total or transport can be present, but if using
# self.correction=="P0", then we should use transport.
if self.correction == "P0" and 'nu-transport' not in self.mgxs_types:
error_flag = True
warn('A "nu-transport" MGXS type is required since a "P0" '
'correction is applied, but a "nu-transport" MGXS is '
'not provided.')
elif self.correction is None and 'total' not in self.mgxs_types:
error_flag = True
warn('A "total" MGXS type is required, but not provided.')
if error_flag:
raise ValueError('Invalid MGXS configuration encountered.')

359
openmc/mgxs/mgxs.py
View file

@ -2717,9 +2717,9 @@ class TransportXS(MGXS):
@property
def scores(self):
if not self.nu:
return ['flux', 'total', 'flux', 'scatter-1']
return ['flux', 'total', 'flux', 'scatter']
else:
return ['flux', 'total', 'flux', 'nu-scatter-1']
return ['flux', 'total', 'flux', 'nu-scatter']
@property
def tally_keys(self):
@ -2730,8 +2730,9 @@ class TransportXS(MGXS):
group_edges = self.energy_groups.group_edges
energy_filter = openmc.EnergyFilter(group_edges)
energyout_filter = openmc.EnergyoutFilter(group_edges)
p1_filter = openmc.LegendreFilter(1)
filters = [[energy_filter], [energy_filter],
[energy_filter], [energyout_filter]]
[energy_filter], [energyout_filter, p1_filter]]
return self._add_angle_filters(filters)
@ -2739,12 +2740,18 @@ class TransportXS(MGXS):
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
# Switch EnergyoutFilter to EnergyFilter.
old_filt = self.tallies['scatter-1'].filters[-1]
p1_tally = self.tallies['scatter-1']
old_filt = p1_tally.filters[-2]
new_filt = openmc.EnergyFilter(old_filt.values)
self.tallies['scatter-1'].filters[-1] = new_filt
p1_tally.filters[-2] = new_filt
self._rxn_rate_tally = \
self.tallies['total'] - self.tallies['scatter-1']
# Slice Legendre expansion filter and change name of score
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
filter_bins=[('P1',)],
squeeze=True)
p1_tally.scores = ['scatter-1']
self._rxn_rate_tally = self.tallies['total'] - p1_tally
self._rxn_rate_tally.sparse = self.sparse
return self._rxn_rate_tally
@ -2758,15 +2765,22 @@ class TransportXS(MGXS):
raise ValueError(msg)
# Switch EnergyoutFilter to EnergyFilter.
old_filt = self.tallies['scatter-1'].filters[-1]
p1_tally = self.tallies['scatter-1']
old_filt = p1_tally.filters[-2]
new_filt = openmc.EnergyFilter(old_filt.values)
self.tallies['scatter-1'].filters[-1] = new_filt
p1_tally.filters[-2] = new_filt
# Slice Legendre expansion filter and change name of score
p1_tally = p1_tally.get_slice(filters=[openmc.LegendreFilter],
filter_bins=[('P1',)],
squeeze=True)
p1_tally.scores = ['scatter-1']
# Compute total cross section
total_xs = self.tallies['total'] / self.tallies['flux (tracklength)']
# Compute transport correction term
trans_corr = self.tallies['scatter-1'] / self.tallies['flux (analog)']
trans_corr = p1_tally / self.tallies['flux (analog)']
# Compute the transport-corrected total cross section
self._xs_tally = total_xs - trans_corr
@ -3510,6 +3524,7 @@ class ScatterXS(MGXS):
self._estimator = 'analog'
self._valid_estimators = ['analog']
class ScatterMatrixXS(MatrixMGXS):
r"""A scattering matrix multi-group cross section with the cosine of the
change-in-angle represented as one or more Legendre moments or a histogram.
@ -3598,10 +3613,10 @@ class ScatterMatrixXS(MatrixMGXS):
name : str, optional
Name of the multi-group cross section. Used as a label to identify
tallies in OpenMC 'tallies.xml' file.
num_polar : Integral, optional
num_polar : int, optional
Number of equi-width polar angle bins for angle discretization;
defaults to one bin
num_azimuthal : Integral, optional
num_azimuthal : int, optional
Number of equi-width azimuthal angle bins for angle discretization;
defaults to one bin
nu : bool
@ -3647,9 +3662,9 @@ class ScatterMatrixXS(MatrixMGXS):
Domain type for spatial homogenization
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
num_polar : Integral
num_polar : int
Number of equi-width polar angle bins for angle discretization
num_azimuthal : Integral
num_azimuthal : int
Number of equi-width azimuthal angle bins for angle discretization
tally_trigger : openmc.Trigger
An (optional) tally precision trigger given to each tally used to
@ -3767,38 +3782,25 @@ class ScatterMatrixXS(MatrixMGXS):
def scores(self):
if self.formulation == 'simple':
scores = ['flux']
if self.scatter_format == 'legendre':
if self.legendre_order == 0:
scores.append('{}-0'.format(self.rxn_type))
if self.correction:
scores.append('{}-1'.format(self.rxn_type))
else:
scores.append('{}-P{}'.format(self.rxn_type, self.legendre_order))
elif self.scatter_format == 'histogram':
scores += [self.rxn_type]
scores = ['flux', self.rxn_type]
else:
# Add scores for groupwise scattering cross section
scores = ['flux', 'scatter']
# Add scores for group-to-group scattering probability matrix
if self.scatter_format == 'legendre':
if self.legendre_order == 0:
scores.append('scatter-0')
else:
scores.append('scatter-P{}'.format(self.legendre_order))
elif self.scatter_format == 'histogram':
scores.append('scatter-0')
# these scores also contain the angular information, whether it be
# Legendre expansion or histogram bins
scores.append('scatter')
# Add scores for multiplicity matrix
# Add scores for multiplicity matrix; scatter info for the
# denominator will come from the previous score
if self.nu:
scores.extend(['nu-scatter-0', 'scatter-0'])
scores.append('nu-scatter')
# Add scores for transport correction
if self.correction == 'P0' and self.legendre_order == 0:
scores.extend(['{}-1'.format(self.rxn_type), 'flux'])
scores.extend([self.rxn_type, 'flux'])
return scores
@ -3811,15 +3813,15 @@ class ScatterMatrixXS(MatrixMGXS):
tally_keys = ['flux (tracklength)', 'scatter']
# Add keys for group-to-group scattering probability matrix
tally_keys.append('scatter-P{}'.format(self.legendre_order))
tally_keys.append('scatter matrix')
# Add keys for multiplicity matrix
if self.nu:
tally_keys.extend(['nu-scatter-0', 'scatter-0'])
tally_keys.extend(['nu-scatter'])
# Add keys for transport correction
if self.correction == 'P0' and self.legendre_order == 0:
tally_keys.extend(['{}-1'.format(self.rxn_type), 'flux (analog)'])
tally_keys.extend(['correction', 'flux (analog)'])
return tally_keys
@ -3836,7 +3838,7 @@ class ScatterMatrixXS(MatrixMGXS):
# Add estimators for multiplicity matrix
if self.nu:
estimators.extend(['analog', 'analog'])
estimators.extend(['analog'])
# Add estimators for transport correction
if self.correction == 'P0' and self.legendre_order == 0:
@ -3853,13 +3855,15 @@ class ScatterMatrixXS(MatrixMGXS):
if self.scatter_format == 'legendre':
if self.correction == 'P0' and self.legendre_order == 0:
filters = [[energy], [energy, energyout], [energyout]]
angle_filter = openmc.LegendreFilter(order=1)
else:
filters = [[energy], [energy, energyout]]
angle_filter = \
openmc.LegendreFilter(order=self.legendre_order)
elif self.scatter_format == 'histogram':
bins = np.linspace(-1., 1., num=self.histogram_bins + 1,
endpoint=True)
filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]]
angle_filter = openmc.MuFilter(bins)
filters = [[energy], [energy, energyout, angle_filter]]
else:
group_edges = self.energy_groups.group_edges
@ -3871,19 +3875,21 @@ class ScatterMatrixXS(MatrixMGXS):
# Group-to-group scattering probability matrix
if self.scatter_format == 'legendre':
filters.append([energy, energyout])
angle_filter = openmc.LegendreFilter(order=self.legendre_order)
elif self.scatter_format == 'histogram':
bins = np.linspace(-1., 1., num=self.histogram_bins + 1,
endpoint=True)
filters.append([energy, energyout, openmc.MuFilter(bins)])
angle_filter = openmc.MuFilter(bins)
filters.append([energy, energyout, angle_filter])
# Multiplicity matrix
if self.nu:
filters.extend([[energy, energyout], [energy, energyout]])
filters.extend([[energy, energyout]])
# Add filters for transport correction
if self.correction == 'P0' and self.legendre_order == 0:
filters.extend([[energyout], [energy]])
filters.extend([[energyout, openmc.LegendreFilter(1)],
[energy]])
return self._add_angle_filters(filters)
@ -3894,27 +3900,39 @@ class ScatterMatrixXS(MatrixMGXS):
if self.formulation == 'simple':
if self.scatter_format == 'legendre':
# If using P0 correction subtract scatter-1 from the diagonal
# If using P0 correction subtract P2 scatter from the diag.
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
scatter_p0 = self.tallies[self.rxn_type].get_slice(
filters=[openmc.LegendreFilter],
filter_bins=[('P0',)])
scatter_p1 = self.tallies[self.rxn_type].get_slice(
filters=[openmc.LegendreFilter],
filter_bins=[('P1',)])
# Transform scatter-p1 tally into an energyin/out matrix
# Set the Legendre order of these tallies to be 0
# so they can be subtracted
legendre = openmc.LegendreFilter(order=0)
scatter_p0.filters[-1] = legendre
scatter_p1.filters[-1] = legendre
scatter_p1 = scatter_p1.summation(
filter_type=openmc.EnergyFilter,
remove_filter=True)
energy_filter = \
scatter_p0.find_filter(openmc.EnergyFilter)
# Transform scatter-p1 into an energyin/out matrix
# to match scattering matrix shape for tally arithmetic
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
scatter_p1 = \
scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = scatter_p0 - scatter_p1
# Extract scattering moment reaction rate Tally
elif self.legendre_order == 0:
tally_key = '{}-{}'.format(self.rxn_type,
self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
# Otherwise, extract scattering moment reaction rate Tally
else:
tally_key = '{}-P{}'.format(self.rxn_type,
self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
self._rxn_rate_tally = self.tallies[self.rxn_type]
elif self.scatter_format == 'histogram':
# Extract scattering rate distribution tally
self._rxn_rate_tally = self.tallies[self.rxn_type]
@ -3941,35 +3959,24 @@ class ScatterMatrixXS(MatrixMGXS):
self._xs_tally = MGXS.xs_tally.fget(self)
else:
# Compute scattering probability matrix
energyout_bins = [self.energy_groups.get_group_bounds(i)
for i in range(self.num_groups, 0, -1)]
tally_key = 'scatter-P{}'.format(self.legendre_order)
# Compute scattering probability matrixS
tally_key = 'scatter matrix'
# Compute normalization factor summed across outgoing energies
norm = self.tallies[tally_key].get_slice(scores=['scatter-0'])
norm = norm.summation(
filter_type=openmc.EnergyoutFilter, filter_bins=energyout_bins)
# Remove the AggregateFilter summed across energyout bins
norm._filters = norm._filters[:2]
if self.scatter_format == 'legendre':
norm = self.tallies[tally_key].get_slice(
scores=['scatter'],
filters=[openmc.LegendreFilter],
filter_bins=[('P0',)], squeeze=True)
# Compute normalization factor summed across outgoing mu bins
if self.scatter_format == 'histogram':
# (Re-)append the MuFilter which was removed above
mu_bins = np.linspace(
-1., 1., num=self.histogram_bins + 1, endpoint=True)
norm._filters.append(openmc.MuFilter(mu_bins))
# Sum across all mu bins
mu_bins = [(mu_bins[i], mu_bins[i+1]) for
i in range(self.histogram_bins)]
elif self.scatter_format == 'histogram':
norm = self.tallies[tally_key].get_slice(
scores=['scatter'])
norm = norm.summation(
filter_type=openmc.MuFilter, filter_bins=mu_bins)
# Remove the AggregateFilter summed across mu bins
norm._filters = norm._filters[:2]
filter_type=openmc.MuFilter, remove_filter=True)
norm = norm.summation(filter_type=openmc.EnergyoutFilter,
remove_filter=True)
# Compute groupwise scattering cross section
self._xs_tally = self.tallies['scatter'] * \
@ -3981,15 +3988,36 @@ class ScatterMatrixXS(MatrixMGXS):
# Multiply by the multiplicity matrix
if self.nu:
numer = self.tallies['nu-scatter-0']
denom = self.tallies['scatter-0']
numer = self.tallies['nu-scatter']
# Get the denominator
if self.scatter_format == 'legendre':
denom = self.tallies[tally_key].get_slice(
scores=['scatter'],
filters=[openmc.LegendreFilter],
filter_bins=[('P0',)], squeeze=True)
# Compute normalization factor summed across mu bins
elif self.scatter_format == 'histogram':
denom = self.tallies[tally_key].get_slice(
scores=['scatter'])
# Sum across all mu bins
denom = denom.summation(
filter_type=openmc.MuFilter, remove_filter=True)
self._xs_tally *= (numer / denom)
# If using P0 correction subtract scatter-1 from the diagonal
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
scatter_p1 = self.tallies['correction'].get_slice(
filters=[openmc.LegendreFilter], filter_bins=[('P1',)])
flux = self.tallies['flux (analog)']
# Set the Legendre order of the P1 tally to be P0
# so it can be subtracted
legendre = openmc.LegendreFilter(order=0)
scatter_p1.filters[-1] = legendre
# Transform scatter-p1 tally into an energyin/out matrix
# to match scattering matrix shape for tally arithmetic
energy_filter = flux.find_filter(openmc.EnergyFilter)
@ -4001,10 +4029,34 @@ class ScatterMatrixXS(MatrixMGXS):
# Override the nuclides for tally arithmetic
correction.nuclides = scatter_p1.nuclides
# Set xs_tally to be itself with only P0 data
self._xs_tally = self._xs_tally.get_slice(
filters=[openmc.LegendreFilter], filter_bins=[('P0',)])
# Tell xs_tally that it is P0
legendre_xs_tally = \
self._xs_tally.find_filter(openmc.LegendreFilter)
legendre_xs_tally.order = 0
# And subtract the P1 correction from the P0 matrix
self._xs_tally -= correction
self._compute_xs()
# Force the angle filter to be the last filter
if self.scatter_format == 'histogram':
angle_filter = self._xs_tally.find_filter(openmc.MuFilter)
else:
angle_filter = \
self._xs_tally.find_filter(openmc.LegendreFilter)
angle_filter_index = self._xs_tally.filters.index(angle_filter)
# If the angle filter index is not last, then make it last
if angle_filter_index != len(self._xs_tally.filters) - 1:
energyout_filter = \
self._xs_tally.find_filter(openmc.EnergyoutFilter)
self._xs_tally._swap_filters(energyout_filter,
angle_filter)
return self._xs_tally
@nu.setter
@ -4125,16 +4177,6 @@ class ScatterMatrixXS(MatrixMGXS):
self._rxn_rate_tally = None
self._loaded_sp = False
if self.scatter_format == 'legendre':
# Expand scores to match the format in the statepoint
# e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
for tally_key, tally in self.tallies.items():
if 'scatter-P' in tally.scores[0]:
score_prefix = tally.scores[0].split('P')[0]
self.tallies[tally_key].scores = \
[score_prefix + '{}'.format(i)
for i in range(self.legendre_order + 1)]
super().load_from_statepoint(statepoint)
def get_slice(self, nuclides=[], in_groups=[], out_groups=[],
@ -4187,12 +4229,11 @@ class ScatterMatrixXS(MatrixMGXS):
slice_xs.legendre_order = legendre_order
# Slice the scattering tally
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
expand_scores = \
[self.rxn_type + '-{}'.format(i)
for i in range(self.legendre_order + 1)]
slice_xs.tallies[tally_key] = \
slice_xs.tallies[tally_key].get_slice(scores=expand_scores)
filter_bins = [tuple(['P{}'.format(i)
for i in range(self.legendre_order + 1)])]
slice_xs.tallies[self.rxn_type] = \
slice_xs.tallies[self.rxn_type].get_slice(
filters=[openmc.LegendreFilter], filter_bins=filter_bins)
# Slice outgoing energy groups if needed
if len(out_groups) != 0:
@ -4206,7 +4247,8 @@ class ScatterMatrixXS(MatrixMGXS):
for tally_type, tally in slice_xs.tallies.items():
if tally.contains_filter(openmc.EnergyoutFilter):
tally_slice = tally.get_slice(
filters=[openmc.EnergyoutFilter], filter_bins=filter_bins)
filters=[openmc.EnergyoutFilter],
filter_bins=filter_bins)
slice_xs.tallies[tally_type] = tally_slice
slice_xs.sparse = self.sparse
@ -4317,14 +4359,19 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins.append((self.energy_groups.get_group_bounds(group),))
# Construct CrossScore for requested scattering moment
if moment != 'all' and self.scatter_format == 'legendre':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
scores = [self.xs_tally.scores[moment]]
if self.scatter_format == 'legendre':
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
filters.append(openmc.LegendreFilter)
filter_bins.append(('P{}'.format(moment),))
num_angle_bins = 1
else:
num_angle_bins = self.legendre_order + 1
else:
scores = []
num_angle_bins = self.histogram_bins
# Construct a collection of the nuclides to retrieve from the xs tally
if self.by_nuclide:
@ -4336,6 +4383,7 @@ class ScatterMatrixXS(MatrixMGXS):
query_nuclides = ['total']
# Use tally summation if user requested the sum for all nuclides
scores = self.xs_tally.scores
if nuclides == 'sum' or nuclides == ['sum']:
xs_tally = self.xs_tally.summation(nuclides=query_nuclides)
xs = xs_tally.get_values(scores=scores, filters=filters,
@ -4367,24 +4415,15 @@ class ScatterMatrixXS(MatrixMGXS):
else:
num_out_groups = len(out_groups)
if self.scatter_format == 'histogram':
num_mu_bins = self.histogram_bins
else:
num_mu_bins = 1
# Reshape tally data array with separate axes for domain and energy
# Accomodate the polar and azimuthal bins if needed
num_subdomains = int(xs.shape[0] / (num_mu_bins * num_in_groups *
num_subdomains = int(xs.shape[0] / (num_angle_bins * num_in_groups *
num_out_groups * self.num_polar *
self.num_azimuthal))
if self.num_polar > 1 or self.num_azimuthal > 1:
if self.scatter_format == 'histogram':
new_shape = (self.num_polar, self.num_azimuthal,
num_subdomains, num_in_groups, num_out_groups,
num_mu_bins)
else:
new_shape = (self.num_polar, self.num_azimuthal,
num_subdomains, num_in_groups, num_out_groups)
new_shape = (self.num_polar, self.num_azimuthal,
num_subdomains, num_in_groups, num_out_groups,
num_angle_bins)
new_shape += xs.shape[1:]
xs = np.reshape(xs, new_shape)
@ -4397,11 +4436,9 @@ class ScatterMatrixXS(MatrixMGXS):
if order_groups == 'increasing':
xs = xs[:, :, :, ::-1, ::-1, ...]
else:
if self.scatter_format == 'histogram':
new_shape = (num_subdomains, num_in_groups, num_out_groups,
num_mu_bins)
else:
new_shape = (num_subdomains, num_in_groups, num_out_groups)
new_shape = (num_subdomains, num_in_groups, num_out_groups,
num_angle_bins)
new_shape += xs.shape[1:]
xs = np.reshape(xs, new_shape)
@ -4416,14 +4453,14 @@ class ScatterMatrixXS(MatrixMGXS):
if squeeze:
# We want to squeeze out everything but the angles, in_groups,
# out_groups, and, if needed, num_mu_bins dimension. These must
# out_groups, and, if needed, num_angle_bins dimension. These must
# not be squeezed so 1-group, 1-angle problems have the correct
# shape.
xs = self._squeeze_xs(xs)
return xs
def get_pandas_dataframe(self, groups='all', nuclides='all', moment='all',
xs_type='macro', paths=True):
def get_pandas_dataframe(self, groups='all', nuclides='all',
xs_type='macro', paths=False):
"""Build a Pandas DataFrame for the MGXS data.
This method leverages :meth:`openmc.Tally.get_pandas_dataframe`, but
@ -4438,19 +4475,15 @@ class ScatterMatrixXS(MatrixMGXS):
may be a list of nuclide name strings (e.g., ['U235', 'U238']).
The special string 'all' will include the cross sections for all
nuclides in the spatial domain. The special string 'sum' will
include the cross sections summed over all nuclides. Defaults
to 'all'.
moment : int or 'all'
The scattering matrix moment to return. All moments will be
returned if the moment is 'all' (default); otherwise, a specific
moment will be returned.
include the cross sections summed over all nuclides. Defaults to
'all'.
xs_type: {'macro', 'micro'}
Return macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
paths : bool, optional
Construct columns for distribcell tally filters (default is True).
The geometric information in the Summary object is embedded into a
Multi-index column with a geometric "path" to each distribcell
The geometric information in the Summary object is embedded into
a Multi-index column with a geometric "path" to each distribcell
instance.
Returns
@ -4466,35 +4499,13 @@ class ScatterMatrixXS(MatrixMGXS):
"""
df = super().get_pandas_dataframe(groups, nuclides, xs_type, paths)
# Build the dataframe using the parent class method
df = super().get_pandas_dataframe(groups, nuclides, xs_type,
paths=paths)
if self.scatter_format == 'legendre':
# Add a moment column to dataframe
if self.legendre_order > 0:
# Insert a column corresponding to the Legendre moments
moments = ['P{}'.format(i)
for i in range(self.legendre_order + 1)]
moments = np.tile(moments, int(df.shape[0] / len(moments)))
df['moment'] = moments
# Place the moment column before the mean column
columns = df.columns.tolist()
mean_index \
= [i for i, s in enumerate(columns) if 'mean' in s][0]
if self.domain_type == 'mesh':
df = df[columns[:mean_index] + [('moment', '')] +
columns[mean_index:-1]]
else:
df = df[columns[:mean_index] + ['moment'] +
columns[mean_index:-1]]
# Select rows corresponding to requested scattering moment
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
df = df[df['moment'] == 'P{}'.format(moment)]
# If the matrix is P0, remove the legendre column
if self.scatter_format == 'legendre' and self.legendre_order == 0:
df = df.drop(axis=1, labels=['legendre'])
return df
@ -4511,8 +4522,9 @@ class ScatterMatrixXS(MatrixMGXS):
The nuclides of the cross-sections to include in the report. This
may be a list of nuclide name strings (e.g., ['U235', 'U238']).
The special string 'all' will report the cross sections for all
nuclides in the spatial domain. The special string 'sum' will report
the cross sections summed over all nuclides. Defaults to 'all'.
nuclides in the spatial domain. The special string 'sum' will
report the cross sections summed over all nuclides. Defaults to
'all'.
xs_type: {'macro', 'micro'}
Return the macro or micro cross section in units of cm^-1 or barns.
Defaults to 'macro'.
@ -4986,14 +4998,9 @@ class ScatterProbabilityMatrix(MatrixMGXS):
def xs_tally(self):
if self._xs_tally is None:
energyout_bins = [self.energy_groups.get_group_bounds(i)
for i in range(self.num_groups, 0, -1)]
norm = self.rxn_rate_tally.get_slice(scores=[self.rxn_type])
norm = norm.summation(
filter_type=openmc.EnergyoutFilter, filter_bins=energyout_bins)
# Remove the AggregateFilter summed across energyout bins
norm._filters = norm._filters[:2]
filter_type=openmc.EnergyoutFilter, remove_filter=True)
# Compute the group-to-group probabilities
self._xs_tally = self.tallies[self.rxn_type] / norm

284
openmc/model/triso.py
View file

@ -113,8 +113,7 @@ class _Domain(metaclass=ABCMeta):
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
Constraint on where particle center can be placed.
volume : float
Volume of the container.
@ -158,8 +157,6 @@ class _Domain(metaclass=ABCMeta):
raise ValueError('Unable to set domain center to {} since it must '
'be of length 3'.format(center))
self._center = [float(x) for x in center]
self._limits = None
self._cell_length = None
def mesh_cell(self, p):
"""Calculate the index of the cell in a mesh overlaid on the domain in
@ -211,6 +208,26 @@ class _Domain(metaclass=ABCMeta):
"""
pass
@abstractmethod
def repel_particles(self, p, q, d, d_new):
"""Move particles p and q apart according to the following
transformation (accounting for boundary conditions on domain):
r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n))
r_j^(n+1) = r_j^(n) - 1/2(d_out^(n+1) - d^(n))
Parameters
----------
p, q : numpy.ndarray
Cartesian coordinates of particle center.
d : float
distance between centers of particles i and j.
d_new : float
final distance between centers of particles i and j.
"""
pass
class _CubicDomain(_Domain):
"""Cubic container in which to pack particles.
@ -238,7 +255,7 @@ class _CubicDomain(_Domain):
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
Maximum distance from center in x-, y-, or z-direction where particle
center can be placed.
volume : float
Volume of the container.
@ -256,9 +273,7 @@ class _CubicDomain(_Domain):
@property
def limits(self):
if self._limits is None:
xlim = self.length/2 - self.particle_radius
self._limits = [[x - xlim for x in self.center],
[x + xlim for x in self.center]]
self._limits = [self.length/2 - self.particle_radius]
return self._limits
@property
@ -284,9 +299,27 @@ class _CubicDomain(_Domain):
self._limits = limits
def random_point(self):
return [uniform(self.limits[0][0], self.limits[1][0]),
uniform(self.limits[0][1], self.limits[1][1]),
uniform(self.limits[0][2], self.limits[1][2])]
x_max = self.limits[0]
return [uniform(-x_max, x_max),
uniform(-x_max, x_max),
uniform(-x_max, x_max)]
def repel_particles(self, p, q, d, d_new):
# Moving each particle distance 's' away from the other along the line
# joining the particle centers will ensure their final distance is
# equal to the outer diameter
s = (d_new - d)/2
v = (p - q)/d
p += s*v
q -= s*v
# Enforce the rigid boundary by moving each particle back along the
# surface normal until it is completely within the container if it
# overlaps the surface
x_max = self.limits[0]
p[:] = np.clip(p, -x_max, x_max)
q[:] = np.clip(q, -x_max, x_max)
class _CylindricalDomain(_Domain):
@ -317,8 +350,8 @@ class _CylindricalDomain(_Domain):
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
Maximum radial distance and maximum distance from center in z-direction
where particle center can be placed.
volume : float
Volume of the container.
@ -340,12 +373,8 @@ class _CylindricalDomain(_Domain):
@property
def limits(self):
if self._limits is None:
xlim = self.length/2 - self.particle_radius
rlim = self.radius - self.particle_radius
self._limits = [[self.center[0] - rlim, self.center[1] - rlim,
self.center[2] - xlim],
[self.center[0] + rlim, self.center[1] + rlim,
self.center[2] + xlim]]
self._limits = [self.radius - self.particle_radius,
self.length/2 - self.particle_radius]
return self._limits
@property
@ -377,10 +406,37 @@ class _CylindricalDomain(_Domain):
self._limits = limits
def random_point(self):
r = sqrt(uniform(0, (self.radius - self.particle_radius)**2))
r_max = self.limits[0]
z_max = self.limits[1]
r = sqrt(uniform(0, r_max**2))
t = uniform(0, 2*pi)
return [r*cos(t) + self.center[0], r*sin(t) + self.center[1],
uniform(self.limits[0][2], self.limits[1][2])]
return [r*cos(t), r*sin(t), uniform(-z_max, z_max)]
def repel_particles(self, p, q, d, d_new):
# Moving each particle distance 's' away from the other along the line
# joining the particle centers will ensure their final distance is
# equal to the outer diameter
s = (d_new - d)/2
v = (p - q)/d
p += s*v
q -= s*v
# Enforce the rigid boundary by moving each particle back along the
# surface normal until it is completely within the container if it
# overlaps the surface
r_max = self.limits[0]
z_max = self.limits[1]
r = sqrt(p[0]**2 + p[1]**2)
if r > r_max:
p[0:2] *= r_max/r
p[2] = np.clip(p[2], -z_max, z_max)
r = sqrt(q[0]**2 + q[1]**2)
if r > r_max:
q[0:2] *= r_max/r
q[2] = np.clip(q[2], -z_max, z_max)
class _SphericalDomain(_Domain):
@ -407,8 +463,7 @@ class _SphericalDomain(_Domain):
Length in x-, y-, and z- directions of each cell in mesh overlaid on
domain.
limits : list of float
Minimum and maximum position in x-, y-, and z-directions where particle
center can be placed.
Maximum radial distance where particle center can be placed.
volume : float
Volume of the container.
@ -425,9 +480,7 @@ class _SphericalDomain(_Domain):
@property
def limits(self):
if self._limits is None:
rlim = self.radius - self.particle_radius
self._limits = [[x - rlim for x in self.center],
[x + rlim for x in self.center]]
self._limits = [self.radius - self.particle_radius]
return self._limits
@property
@ -453,10 +506,33 @@ class _SphericalDomain(_Domain):
self._limits = limits
def random_point(self):
r_max = self.limits[0]
x = (gauss(0, 1), gauss(0, 1), gauss(0, 1))
r = (uniform(0, (self.radius - self.particle_radius)**3)**(1/3) /
sqrt(x[0]**2 + x[1]**2 + x[2]**2))
return [r*x[i] + self.center[i] for i in range(3)]
r = (uniform(0, r_max**3)**(1/3) / sqrt(x[0]**2 + x[1]**2 + x[2]**2))
return [r*s for s in x]
def repel_particles(self, p, q, d, d_new):
# Moving each particle distance 's' away from the other along the line
# joining the particle centers will ensure their final distance is
# equal to the outer diameter
s = (d_new - d)/2
v = (p - q)/d
p += s*v
q -= s*v
# Enforce the rigid boundary by moving each particle back along the
# surface normal until it is completely within the container if it
# overlaps the surface
r_max = self.limits[0]
r = sqrt(p[0]**2 + p[1]**2 + p[2]**2)
if r > r_max:
p *= r_max/r
r = sqrt(q[0]**2 + q[1]**2 + q[2]**2)
if r > r_max:
q *= r_max/r
def create_triso_lattice(trisos, lower_left, pitch, shape, background):
@ -636,6 +712,7 @@ def _close_random_pack(domain, particles, contraction_rate):
del rods_map[i]
del rods_map[j]
return d, i, j
return None, None, None
def create_rod_list():
"""Generate sorted list of rods (distances between particle centers).
@ -660,8 +737,8 @@ def _close_random_pack(domain, particles, contraction_rate):
# Find distance to nearest neighbor and index of nearest neighbor for
# all particles
d, n = tree.query(particles, k=2)
d = d[:,1]
n = n[:,1]
d = d[:, 1]
n = n[:, 1]
# Array of particle indices, indices of nearest neighbors, and
# distances to nearest neighbors
@ -670,8 +747,8 @@ def _close_random_pack(domain, particles, contraction_rate):
# Sort along second column and swap first and second columns to create
# array of nearest neighbor indices, indices of particles they are
# nearest neighbors of, and distances between them
b = a[a[:,1].argsort()]
b[:,[0, 1]] = b[:,[1, 0]]
b = a[a[:, 1].argsort()]
b[:, [0, 1]] = b[:, [1, 0]]
# Find the intersection between 'a' and 'b': a list of particles who
# are each other's nearest neighbors and the distance between them
@ -685,12 +762,8 @@ def _close_random_pack(domain, particles, contraction_rate):
del rods[:]
rods_map.clear()
for d, i, j in r:
add_rod(d, i, j)
# Inner diameter is set initially to the shortest center-to-center
# distance between any two particles
if rods:
inner_diameter[0] = rods[0][0]
if d < outer_diameter and not np.isclose(d, outer_diameter, atol=1.0e-14):
add_rod(d, i, j)
def update_mesh(i):
"""Update which mesh cells the particle is in based on new particle
@ -729,50 +802,19 @@ def _close_random_pack(domain, particles, contraction_rate):
j = floor(-log10(pf_out - pf_in)).
Returns
-------
float
New outer diameter
"""
inner_pf = (4/3 * pi * (inner_diameter[0]/2)**3 * n_particles /
domain.volume)
outer_pf = (4/3 * pi * (outer_diameter[0]/2)**3 * n_particles /
domain.volume)
inner_pf = 4/3*pi*(inner_diameter/2)**3*n_particles/domain.volume
outer_pf = 4/3*pi*(outer_diameter/2)**3*n_particles/domain.volume
j = floor(-log10(outer_pf - inner_pf))
outer_diameter[0] = (outer_diameter[0] - 0.5**j * contraction_rate *
initial_outer_diameter / n_particles)
def repel_particles(i, j, d):
"""Move particles p and q apart according to the following
transformation (accounting for reflective boundary conditions on
domain):
r_i^(n+1) = r_i^(n) + 1/2(d_out^(n+1) - d^(n))
r_j^(n+1) = r_j^(n) - 1/2(d_out^(n+1) - d^(n))
Parameters
----------
i, j : int
Index of particles in particles array.
d : float
distance between centers of particles i and j.
"""
# Moving each particle distance 'r' away from the other along the line
# joining the particle centers will ensure their final distance is equal
# to the outer diameter
r = (outer_diameter[0] - d)/2
v = (particles[i] - particles[j])/d
particles[i] += r*v
particles[j] -= r*v
# Apply reflective boundary conditions
particles[i] = particles[i].clip(domain.limits[0], domain.limits[1])
particles[j] = particles[j].clip(domain.limits[0], domain.limits[1])
update_mesh(i)
update_mesh(j)
return (outer_diameter - 0.5**j * contraction_rate *
initial_outer_diameter / n_particles)
def nearest(i):
"""Find index of nearest neighbor of particle i.
@ -803,14 +845,14 @@ def _close_random_pack(domain, particles, contraction_rate):
else:
return None, None
def update_rod_list(i, j):
"""Update the rod list with the new nearest neighbors of particles i
and j since their overlap was eliminated.
def update_rod_list(i):
"""Update the rod list with the new nearest neighbors of particle since
its overlap was eliminated.
Parameters
----------
i, j : int
Index of particles in particles array.
i : int
Index of particle in particles array.
"""
@ -818,18 +860,10 @@ def _close_random_pack(domain, particles, contraction_rate):
# remove the rod currently containing k from the rod list and add rod
# k-i, keeping the rod list sorted
k, d_ik = nearest(i)
if k and nearest(k)[0] == i:
if (k and nearest(k)[0] == i and d_ik < outer_diameter
and not np.isclose(d, outer_diameter, atol=1.0e-14)):
remove_rod(k)
add_rod(d_ik, i, k)
l, d_jl = nearest(j)
if l and nearest(l)[0] == j:
remove_rod(l)
add_rod(d_jl, j, l)
# Set inner diameter to the shortest distance between two particle
# centers
if rods:
inner_diameter[0] = rods[0][0]
n_particles = len(particles)
diameter = 2*domain.particle_radius
@ -841,36 +875,71 @@ def _close_random_pack(domain, particles, contraction_rate):
initial_outer_diameter = 2*(domain.volume/(n_particles*4/3*pi))**(1/3)
# Inner and outer diameter of particles will change during packing
outer_diameter = [initial_outer_diameter]
inner_diameter = [0]
outer_diameter = initial_outer_diameter
inner_diameter = 0.
# List of rods arranged in a heap and mapping of particle ids to rods
rods = []
rods_map = {}
# Initialize two-way dictionary that identifies which particles are near a
# given mesh cell and which mesh cells a particle is near
mesh = defaultdict(set)
mesh_map = defaultdict(set)
for i in range(n_particles):
for idx in domain.nearby_mesh_cells(particles[i]):
mesh[idx].add(i)
mesh_map[i].add(idx)
while True:
# Rebuild the sorted list of rods according to the current particle
# configuration
create_rod_list()
if inner_diameter[0] >= diameter:
# Set the inner diameter to the shortest center-to-center distance
# between any two particles
if rods:
inner_diameter = rods[0][0]
# Reached the desired particle radius
if inner_diameter >= diameter:
break
# The algorithm converged before reaching the desired particle radius.
# This can happen when the desired packing fraction is close to the
# packing fraction limit. The packing fraction is a random variable
# that is determined by the particle locations and the contraction
# rate. A higher packing fraction can be achieved with a smaller
# contraction rate, though at the cost of a longer simulation time --
# the number of iterations needed to remove all overlaps is inversely
# proportional to the contraction rate.
if inner_diameter >= outer_diameter or not rods:
warnings.warn('Close random pack converged before reaching true '
'particle radius; some particles may overlap. Try '
'reducing contraction rate or packing fraction.')
break
while True:
d, i, j = pop_rod()
reduce_outer_diameter()
repel_particles(i, j, d)
update_rod_list(i, j)
if inner_diameter[0] >= diameter or not rods:
if not d:
break
outer_diameter = reduce_outer_diameter()
domain.repel_particles(particles[i], particles[j], d, outer_diameter)
update_mesh(i)
update_mesh(j)
update_rod_list(i)
update_rod_list(j)
if not rods:
break
inner_diameter = rods[0][0]
if inner_diameter >= diameter or inner_diameter >= outer_diameter:
break
def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
domain_radius=None, domain_center=[0., 0., 0.],
n_particles=None, packing_fraction=None,
initial_packing_fraction=0.3, contraction_rate=1/400, seed=1):
initial_packing_fraction=0.3, contraction_rate=1.e-3, seed=1):
"""Generate a random, non-overlapping configuration of TRISO particles
within a container.
@ -933,7 +1002,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
to speed up the nearest neighbor search by only searching for a particle's
neighbors within that mesh cell.
In CRP, each particle is assigned two diameters, and inner and an outer,
In CRP, each particle is assigned two diameters, an inner and an outer,
which approach each other during the simulation. The inner diameter,
defined as the minimum center-to-center distance, is the true diameter of
the particles and defines the pf. At each iteration the worst overlap
@ -1008,8 +1077,7 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
# Recalculate the limits for the initial random sequential packing using
# the desired final particle radius to ensure particles are fully contained
# within the domain during the close random pack
domain.limits = [[x - initial_radius + radius for x in domain.limits[0]],
[x + initial_radius - radius for x in domain.limits[1]]]
domain.limits = [x + initial_radius - radius for x in domain.limits]
# Generate non-overlapping particles for an initial inner radius using
# random sequential packing algorithm
@ -1024,5 +1092,5 @@ def pack_trisos(radius, fill, domain_shape='cylinder', domain_length=None,
trisos = []
for p in particles:
trisos.append(TRISO(radius, fill, p))
trisos.append(TRISO(radius, fill, [x + c for x, c in zip(p, domain.center)]))
return trisos

23
openmc/plotter.py
View file

@ -170,13 +170,13 @@ def plot_xs(this, types, divisor_types=None, temperature=294., data_type=None,
data = data_new
else:
# Calculate for MG cross sections
E, data = calculate_mgxs(this, types, orders, temperature,
E, data = calculate_mgxs(this, data_type, types, orders, temperature,
mg_cross_sections, ce_cross_sections,
enrichment)
if divisor_types:
cv.check_length('divisor types', divisor_types, len(types))
Ediv, data_div = calculate_mgxs(this, divisor_types,
Ediv, data_div = calculate_mgxs(this, data_type, divisor_types,
divisor_orders, temperature,
mg_cross_sections,
ce_cross_sections, enrichment)
@ -243,7 +243,7 @@ def calculate_cexs(this, data_type, types, temperature=294., sab_name=None,
Parameters
----------
this : str or openmc.Material
this : {str, openmc.Nuclide, openmc.Element, openmc.Material}
Object to source data from
data_type : {'nuclide', 'element', material'}
Type of object to plot
@ -280,7 +280,11 @@ def calculate_cexs(this, data_type, types, temperature=294., sab_name=None,
cv.check_type('enrichment', enrichment, Real)
if data_type == 'nuclide':
energy_grid, xs = _calculate_cexs_nuclide(this, types, temperature,
if isinstance(this, str):
nuc = openmc.Nuclide(this)
else:
nuc = this
energy_grid, xs = _calculate_cexs_nuclide(nuc, types, temperature,
sab_name, cross_sections)
# Convert xs (Iterable of Callable) to a grid of cross section values
# calculated on @ the points in energy_grid for consistency with the
@ -289,10 +293,15 @@ def calculate_cexs(this, data_type, types, temperature=294., sab_name=None,
for line in range(len(types)):
data[line, :] = xs[line](energy_grid)
elif data_type == 'element':
energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,
if isinstance(this, str):
elem = openmc.Element(this)
else:
elem = this
energy_grid, data = _calculate_cexs_elem_mat(elem, types, temperature,
cross_sections, sab_name,
enrichment)
elif data_type == 'material':
cv.check_type('this', this, openmc.Material)
energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,
cross_sections)
else:
@ -518,10 +527,8 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
T = this.temperature
else:
T = temperature
data_type = 'material'
else:
T = temperature
data_type = 'element'
# Load the library
library = openmc.data.DataLibrary.from_xml(cross_sections)
@ -571,7 +578,7 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
name = nuclide[0]
nuc = nuclide[1]
sab_tab = sabs[name]
temp_E, temp_xs = calculate_cexs(nuc, data_type, types, T, sab_tab,
temp_E, temp_xs = calculate_cexs(nuc, 'nuclide', types, T, sab_tab,
cross_sections)
E.append(temp_E)
# Since the energy grids are different, store the cross sections as

9
openmc/statepoint.py
View file

@ -150,6 +150,8 @@ class StatePoint(object):
def __exit__(self, *exc):
self._f.close()
if self._summary is not None:
self._summary._f.close()
@property
def cmfd_on(self):
@ -405,17 +407,10 @@ class StatePoint(object):
scores = group['score_bins'].value
n_score_bins = group['n_score_bins'].value
# Read scattering moment order strings (e.g., P3, Y1,2, etc.)
moments = group['moment_orders'].value
# Add the scores to the Tally
for j, score in enumerate(scores):
score = score.decode()
# If this is a moment, use generic moment order
pattern = r'-n$|-pn$|-yn$'
score = re.sub(pattern, '-' + moments[j].decode(), score)
tally.scores.append(score)
# Add Tally to the global dictionary of all Tallies

27
openmc/summary.py
View file

@ -9,7 +9,7 @@ import openmc
import openmc.checkvalue as cv
from openmc.region import Region
_VERSION_SUMMARY = 5
_VERSION_SUMMARY = 6
class Summary(object):
@ -26,6 +26,8 @@ class Summary(object):
nuclides : dict
Dictionary whose keys are nuclide names and values are atomic weight
ratios.
macroscopics : list
Names of macroscopic data sets
version: tuple of int
Version of OpenMC
@ -44,13 +46,15 @@ class Summary(object):
self._fast_materials = {}
self._fast_surfaces = {}
self._fast_cells = {}
self._fast_universes = {}
self._fast_universes = {}
self._fast_lattices = {}
self._materials = openmc.Materials()
self._nuclides = {}
self._macroscopics = []
self._read_nuclides()
self._read_macroscopics()
with warnings.catch_warnings():
warnings.simplefilter("ignore", openmc.IDWarning)
self._read_geometry()
@ -71,15 +75,26 @@ class Summary(object):
def nuclides(self):
return self._nuclides
@property
def macroscopics(self):
return self._macroscopics
@property
def version(self):
return tuple(self._f.attrs['openmc_version'])
def _read_nuclides(self):
names = self._f['nuclides/names'].value
awrs = self._f['nuclides/awrs'].value
for name, awr in zip(names, awrs):
self._nuclides[name.decode()] = awr
if 'nuclides/names' in self._f:
names = self._f['nuclides/names'].value
awrs = self._f['nuclides/awrs'].value
for name, awr in zip(names, awrs):
self._nuclides[name.decode()] = awr
def _read_macroscopics(self):
if 'macroscopics/names' in self._f:
names = self._f['macroscopics/names'].value
for name in names:
self._macroscopics = name.decode()
def _read_geometry(self):
# Read in and initialize the Materials and Geometry

75
openmc/tallies.py
View file

@ -65,9 +65,7 @@ class Tally(IDManagerMixin):
triggers : list of openmc.Trigger
List of tally triggers
num_scores : int
Total number of scores, accounting for the fact that a single
user-specified score, e.g. scatter-P3 or flux-Y2,2, might have multiple
bins
Total number of scores
num_filter_bins : int
Total number of filter bins accounting for all filters
num_bins : int
@ -388,6 +386,13 @@ class Tally(IDManagerMixin):
# If score is a string, strip whitespace
if isinstance(score, str):
# Check to see if scores are deprecated before storing
for deprecated in ['scatter-', 'nu-scatter-', 'scatter-p',
'nu-scatter-p', 'scatter-y', 'nu-scatter-y',
'flux-y', 'total-y']:
if score.startswith(deprecated):
msg = score.strip() + ' is no longer supported.'
raise ValueError(msg)
scores[i] = score.strip()
self._scores = cv.CheckedList(_SCORE_CLASSES, 'tally scores', scores)
@ -827,51 +832,8 @@ class Tally(IDManagerMixin):
# Sparsify merged tally if both tallies are sparse
merged_tally.sparse = self.sparse and other.sparse
# Consolidate scatter and flux Legendre moment scores
merged_tally._consolidate_moment_scores()
return merged_tally
def _consolidate_moment_scores(self):
"""Remove redundant scattering and flux moment scores from a Tally."""
# Define regex for scatter, nu-scatter and flux moment scores
regex = [(r'^((?!nu-)scatter-\d)', r'^((?!nu-)scatter-(P|p)\d)'),
(r'nu-scatter-\d', r'nu-scatter-(P|p)\d'),
(r'flux-\d', r'flux-(P|p)\d')]
# Find all non-scattering and non-flux moment scores
scores = [x for x in self.scores if
re.search(r'^((?!scatter-).)*$', x)]
scores = [x for x in scores if
re.search(r'^((?!flux-).)*$', x)]
for regex_n, regex_pn in regex:
# Use regex to find score-(P)n scores
score_n = [x for x in self.scores if re.search(regex_n, x)]
score_pn = [x for x in self.scores if re.search(regex_pn, x)]
# Consolidate moment scores
if len(score_pn) > 0:
# Only keep the highest score-PN score
high_pn = sorted([x.lower() for x in score_pn])[-1]
pn = int(high_pn.split('-')[-1].replace('p', ''))
# Only keep the score-N scores with N > PN
score_n = sorted([x.lower() for x in score_n])
score_n = [x for x in score_n if (int(x.split('-')[1]) > pn)]
# Append highest score-PN and any higher score-N scores
scores.extend([high_pn] + score_n)
else:
scores.extend(score_n)
# Override Tally's scores with consolidated list of scores
self.scores = scores
def to_xml_element(self):
"""Return XML representation of the tally
@ -1180,7 +1142,7 @@ class Tally(IDManagerMixin):
# Determine the score indices from any of the requested scores
if nuclides:
nuclide_indices = np.zeros(len(nuclides), dtype=np.int)
nuclide_indices = np.zeros(len(nuclides), dtype=int)
for i, nuclide in enumerate(nuclides):
nuclide_indices[i] = self.get_nuclide_index(nuclide)
@ -1219,7 +1181,7 @@ class Tally(IDManagerMixin):
# Determine the score indices from any of the requested scores
if scores:
score_indices = np.zeros(len(scores), dtype=np.int)
score_indices = np.zeros(len(scores), dtype=int)
for i, score in enumerate(scores):
score_indices[i] = self.get_score_index(score)
@ -1491,11 +1453,8 @@ class Tally(IDManagerMixin):
data = self.get_values(value=value)
# Build a new array shape with one dimension per filter
new_shape = ()
for self_filter in self.filters:
new_shape += (self_filter.num_bins, )
new_shape += (self.num_nuclides,)
new_shape += (self.num_scores,)
new_shape = tuple(f.num_bins for f in self.filters)
new_shape += (self.num_nuclides, self.num_scores)
# Reshape the data with one dimension for each filter
data = np.reshape(data, new_shape)
@ -2772,7 +2731,7 @@ class Tally(IDManagerMixin):
# Sum across the bins in the user-specified filter
for i, self_filter in enumerate(self.filters):
if isinstance(self_filter, filter_type):
if type(self_filter) == filter_type:
shape = mean.shape
mean = np.take(mean, indices=bin_indices, axis=i)
std_dev = np.take(std_dev, indices=bin_indices, axis=i)
@ -3012,7 +2971,7 @@ class Tally(IDManagerMixin):
The data in the derived tally arrays is "diagonalized" along the bins in
the new filter. This functionality is used by the openmc.mgxs module; to
transport-correct scattering matrices by subtracting a 'scatter-P1'
reaction rate tally with an energy filter from an 'scatter' reaction
reaction rate tally with an energy filter from a 'scatter' reaction
rate tally with both energy and energyout filters.
Parameters
@ -3031,7 +2990,7 @@ class Tally(IDManagerMixin):
if new_filter in self.filters:
msg = 'Unable to diagonalize Tally ID="{0}" which already ' \
'contains a "{1}" filter'.format(self.id, new_filter.type)
'contains a "{1}" filter'.format(self.id, type(new_filter))
raise ValueError(msg)
# Add the new filter to a copy of this Tally
@ -3042,8 +3001,8 @@ class Tally(IDManagerMixin):
# by which the "base" indices should be repeated to account for all
# other filter bins in the diagonalized tally
indices = np.arange(0, new_filter.num_bins**2, new_filter.num_bins+1)
diag_factor = int(self.num_filter_bins / new_filter.num_bins)
diag_indices = np.zeros(self.num_filter_bins, dtype=np.int)
diag_factor = self.num_filter_bins // new_filter.num_bins
diag_indices = np.zeros(self.num_filter_bins, dtype=int)
# Determine the filter indices along the new "diagonal"
for i in range(diag_factor):

4
src/angle_distribution.F90
View file

@ -1,12 +1,10 @@
module angle_distribution
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: read_attribute, get_shape, read_dataset, &
open_dataset, close_dataset
open_dataset, close_dataset, HID_T, HSIZE_T
use random_lcg, only: prn
implicit none

2
src/angleenergy_header.F90
View file

@ -1,6 +1,6 @@
module angleenergy_header
use hdf5, only: HID_T
use hdf5_interface, only: HID_T
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy

42
src/api.F90
View file

@ -2,8 +2,6 @@ module openmc_api
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, h5tclose_f, h5close_f
use bank_header, only: openmc_source_bank
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum, openmc_get_keff
@ -15,7 +13,7 @@ module openmc_api
use mesh_header
use message_passing
use nuclide_header
use initialize, only: openmc_init
use initialize, only: openmc_init_f
use particle_header, only: Particle
use plot, only: openmc_plot_geometry
use random_lcg, only: openmc_get_seed, openmc_set_seed
@ -64,7 +62,7 @@ module openmc_api
public :: openmc_get_tally_index
public :: openmc_global_tallies
public :: openmc_hard_reset
public :: openmc_init
public :: openmc_init_f
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
@ -80,7 +78,6 @@ module openmc_api
public :: openmc_nuclide_name
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_run
public :: openmc_set_seed
public :: openmc_simulation_finalize
public :: openmc_simulation_init
@ -105,12 +102,16 @@ contains
! variables
!===============================================================================
subroutine openmc_finalize() bind(C)
function openmc_finalize() result(err) bind(C)
integer(C_INT) :: err
integer :: err
interface
subroutine openmc_free_bank() bind(C)
end subroutine openmc_free_bank
end interface
! Clear results
call openmc_reset()
err = openmc_reset()
! Reset global variables
assume_separate = .false.
@ -170,18 +171,14 @@ contains
! Deallocate arrays
call free_memory()
! Release compound datatypes
call h5tclose_f(hdf5_bank_t, err)
! Close FORTRAN interface.
call h5close_f(err)
err = 0
#ifdef OPENMC_MPI
! Free all MPI types
call MPI_TYPE_FREE(MPI_BANK, err)
call openmc_free_bank()
#endif
end subroutine openmc_finalize
end function openmc_finalize
!===============================================================================
! OPENMC_FIND determines the ID or a cell or material at a given point in space
@ -232,9 +229,11 @@ contains
! generator state
!===============================================================================
subroutine openmc_hard_reset() bind(C)
function openmc_hard_reset() result(err) bind(C)
integer(C_INT) :: err
! Reset all tallies and timers
call openmc_reset()
err = openmc_reset()
! Reset total generations and keff guess
keff = ONE
@ -242,13 +241,15 @@ contains
! Reset the random number generator state
call openmc_set_seed(DEFAULT_SEED)
end subroutine openmc_hard_reset
end function openmc_hard_reset
!===============================================================================
! OPENMC_RESET resets tallies and timers
!===============================================================================
subroutine openmc_reset() bind(C)
function openmc_reset() result(err) bind(C)
integer(C_INT) :: err
integer :: i
if (allocated(tallies)) then
@ -296,7 +297,8 @@ contains
call time_transport % reset()
call time_finalize % reset()
end subroutine openmc_reset
err = 0
end function openmc_reset
!===============================================================================
! FREE_MEMORY deallocates and clears all global allocatable arrays in the

53
src/cmfd_data.F90
View file

@ -76,6 +76,7 @@ contains
integer :: i_filter_mesh ! index for mesh filter
integer :: i_filter_ein ! index for incoming energy filter
integer :: i_filter_eout ! index for outgoing energy filter
integer :: i_filter_legendre ! index for Legendre filter
integer :: i_mesh ! flattend index for mesh
logical :: energy_filters! energy filters present
real(8) :: flux ! temp variable for flux
@ -116,8 +117,11 @@ contains
if (ital < 3) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
else
else if (ital == 3) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESHSURFACE))
else if (ital == 4) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_legendre = t % filter(t % find_filter(FILTER_LEGENDRE))
end if
! Check for energy filters
@ -187,13 +191,6 @@ contains
! Get total rr and convert to total xs
cmfd % totalxs(h,i,j,k) = t % results(RESULT_SUM,2,score_index) / flux
! Get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = t % results(RESULT_SUM,3,score_index) / flux
! Calculate diffusion coefficient
cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - &
cmfd % p1scattxs(h,i,j,k)))
else if (ital == 2) then
! Begin loop to get energy out tallies
@ -301,6 +298,46 @@ contains
score_index + IN_TOP)
cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + OUT_TOP)
else if (ital == 4) then
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins % clear()
call filter_matches(t % filter(l)) % bins % push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
filter_matches(i_filter_mesh) % bins % data(1) = &
m % get_bin_from_indices(ijk)
! Apply energy in filter
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
end if
! Apply Legendre filter
filter_matches(i_filter_legendre) % bins % data(1) = 2
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins % data(1) - 1) * t % stride(l)
end do
! Get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = &
t % results(RESULT_SUM,1,score_index) / &
cmfd % flux(h,i,j,k)
! Calculate diffusion coefficient
cmfd % diffcof(h,i,j,k) = &
ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - &
cmfd % p1scattxs(h,i,j,k)))
end if TALLY
end do OUTGROUP

68
src/cmfd_input.F90
View file

@ -278,7 +278,7 @@ contains
m % id = i_start
! Set mesh type to rectangular
m % type = LATTICE_RECT
m % type = MESH_REGULAR
! Get pointer to mesh XML node
node_mesh = root % child("mesh")
@ -373,7 +373,7 @@ contains
! Determine number of filters
energy_filters = check_for_node(node_mesh, "energy")
n = merge(4, 2, energy_filters)
n = merge(5, 3, energy_filters)
! Extend filters array so we can add CMFD filters
err = openmc_extend_filters(n, i_filt_start, i_filt_end)
@ -414,6 +414,12 @@ contains
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_meshsurface_filter_set_mesh(i_filt, i_start)
! Add in legendre filter for the P1 tally
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'legendre' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_legendre_filter_set_order(i_filt, 1)
! Initialize filters
do i = i_filt_start, i_filt_end
@ -421,7 +427,7 @@ contains
end do
! Allocate tallies
err = openmc_extend_tallies(3, i_start, i_end)
err = openmc_extend_tallies(4, i_start, i_end)
cmfd_tallies => tallies(i_start:i_end)
! Begin loop around tallies
@ -455,7 +461,7 @@ contains
if (i == 1) then
! Set name
t % name = "CMFD flux, total, scatter-1"
t % name = "CMFD flux, total"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
@ -473,19 +479,12 @@ contains
deallocate(filter_indices)
! Allocate scoring bins
allocate(t % score_bins(3))
t % n_score_bins = 3
t % n_user_score_bins = 3
! Allocate scattering order data
allocate(t % moment_order(3))
t % moment_order = 0
allocate(t % score_bins(2))
t % n_score_bins = 2
! Set macro_bins
t % score_bins(1) = SCORE_FLUX
t % score_bins(2) = SCORE_TOTAL
t % score_bins(3) = SCORE_SCATTER_N
t % moment_order(3) = 1
else if (i == 2) then
@ -517,11 +516,6 @@ contains
! Allocate macro reactions
allocate(t % score_bins(2))
t % n_score_bins = 2
t % n_user_score_bins = 2
! Allocate scattering order data
allocate(t % moment_order(2))
t % moment_order = 0
! Set macro_bins
t % score_bins(1) = SCORE_NU_SCATTER
@ -537,7 +531,7 @@ contains
! Allocate and set filters
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_end
filter_indices(1) = i_filt_end - 1
if (energy_filters) then
filter_indices(2) = i_filt_start + 1
end if
@ -547,15 +541,41 @@ contains
! Allocate macro reactions
allocate(t % score_bins(1))
t % n_score_bins = 1
t % n_user_score_bins = 1
! Allocate scattering order data
allocate(t % moment_order(1))
t % moment_order = 0
! Set macro bins
t % score_bins(1) = SCORE_CURRENT
t % type = TALLY_MESH_SURFACE
else if (i == 4) then
! Set name
t % name = "CMFD P1 scatter"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
! Set tally type to volume
t % type = TALLY_VOLUME
! Allocate and set filters
n_filter = 2
if (energy_filters) then
n_filter = n_filter + 1
end if
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_start
filter_indices(2) = i_filt_end
if (energy_filters) then
filter_indices(3) = i_filt_start + 1
end if
err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices)
deallocate(filter_indices)
! Allocate scoring bins
allocate(t % score_bins(1))
t % n_score_bins = 1
! Set macro_bins
t % score_bins(1) = SCORE_SCATTER
end if
! Make CMFD tallies active from the start

54
src/constants.F90
View file

@ -21,7 +21,7 @@ module constants
integer, parameter :: VERSION_STATEPOINT(2) = [17, 0]
integer, parameter :: VERSION_PARTICLE_RESTART(2) = [2, 0]
integer, parameter :: VERSION_TRACK(2) = [2, 0]
integer, parameter :: VERSION_SUMMARY(2) = [5, 0]
integer, parameter :: VERSION_SUMMARY(2) = [6, 0]
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0]
@ -232,6 +232,10 @@ module constants
MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data
MGXS_ANGLE = 2 ! Data by Angular Bins
! Flag to denote this was a macroscopic data object
real(8), parameter :: &
MACROSCOPIC_AWR = -TWO
! Fission neutron emission (nu) type
integer, parameter :: &
NU_NONE = 0, & ! No nu values (non-fissionable)
@ -310,50 +314,28 @@ module constants
! Tally score type -- if you change these, make sure you also update the
! _SCORES dictionary in openmc/capi/tally.py
integer, parameter :: N_SCORE_TYPES = 24
integer, parameter :: N_SCORE_TYPES = 16
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
SCORE_SCATTER = -3, & ! scattering rate
SCORE_NU_SCATTER = -4, & ! scattering production rate
SCORE_SCATTER_N = -5, & ! arbitrary scattering moment
SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment
SCORE_NU_SCATTER_N = -7, & ! arbitrary nu-scattering moment
SCORE_NU_SCATTER_PN = -8, & ! system for scoring 0th through nth nu-scatter moment
SCORE_ABSORPTION = -9, & ! absorption rate
SCORE_FISSION = -10, & ! fission rate
SCORE_NU_FISSION = -11, & ! neutron production rate
SCORE_KAPPA_FISSION = -12, & ! fission energy production rate
SCORE_CURRENT = -13, & ! current
SCORE_FLUX_YN = -14, & ! angular moment of flux
SCORE_TOTAL_YN = -15, & ! angular moment of total reaction rate
SCORE_SCATTER_YN = -16, & ! angular flux-weighted scattering moment (0:N)
SCORE_NU_SCATTER_YN = -17, & ! angular flux-weighted nu-scattering moment (0:N)
SCORE_EVENTS = -18, & ! number of events
SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate
SCORE_PROMPT_NU_FISSION = -20, & ! prompt neutron production rate
SCORE_INVERSE_VELOCITY = -21, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value
SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate
SCORE_ABSORPTION = -5, & ! absorption rate
SCORE_FISSION = -6, & ! fission rate
SCORE_NU_FISSION = -7, & ! neutron production rate
SCORE_KAPPA_FISSION = -8, & ! fission energy production rate
SCORE_CURRENT = -9, & ! current
SCORE_EVENTS = -10, & ! number of events
SCORE_DELAYED_NU_FISSION = -11, & ! delayed neutron production rate
SCORE_PROMPT_NU_FISSION = -12, & ! prompt neutron production rate
SCORE_INVERSE_VELOCITY = -13, & ! flux-weighted inverse velocity
SCORE_FISS_Q_PROMPT = -14, & ! prompt fission Q-value
SCORE_FISS_Q_RECOV = -15, & ! recoverable fission Q-value
SCORE_DECAY_RATE = -16 ! delayed neutron precursor decay rate
! Maximum scattering order supported
integer, parameter :: MAX_ANG_ORDER = 10
! Names of *-PN & *-YN scores (MOMENT_STRS) and *-N moment scores
character(*), parameter :: &
MOMENT_STRS(6) = (/ "scatter-p ", &
"nu-scatter-p", &
"flux-y ", &
"total-y ", &
"scatter-y ", &
"nu-scatter-y"/), &
MOMENT_N_STRS(2) = (/ "scatter- ", &
"nu-scatter- "/)
! Location in MOMENT_STRS where the YN data begins
integer, parameter :: YN_LOC = 3
! Tally map bin finding
integer, parameter :: NO_BIN_FOUND = -1

16
src/endf.F90
View file

@ -26,14 +26,6 @@ contains
string = "scatter"
case (SCORE_NU_SCATTER)
string = "nu-scatter"
case (SCORE_SCATTER_N)
string = "scatter-n"
case (SCORE_SCATTER_PN)
string = "scatter-pn"
case (SCORE_NU_SCATTER_N)
string = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
string = "nu-scatter-pn"
case (SCORE_ABSORPTION)
string = "absorption"
case (SCORE_FISSION)
@ -50,14 +42,6 @@ contains
string = "kappa-fission"
case (SCORE_CURRENT)
string = "current"
case (SCORE_FLUX_YN)
string = "flux-yn"
case (SCORE_TOTAL_YN)
string = "total-yn"
case (SCORE_SCATTER_YN)
string = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
string = "nu-scatter-yn"
case (SCORE_EVENTS)
string = "events"
case (SCORE_INVERSE_VELOCITY)

2
src/endf_header.F90
View file

@ -1,7 +1,5 @@
module endf_header
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use constants, only: ZERO, HISTOGRAM, LINEAR_LINEAR, LINEAR_LOG, &
LOG_LINEAR, LOG_LOG

2
src/energy_distribution.F90
View file

@ -1,7 +1,5 @@
module energy_distribution
use hdf5
use algorithm, only: binary_search
use constants, only: ZERO, ONE, HALF, TWO, PI, HISTOGRAM, LINEAR_LINEAR
use endf_header, only: Tabulated1D

30
src/error.F90
View file

@ -15,19 +15,19 @@ module error
public :: write_message
! Error codes
integer(C_INT), public, bind(C) :: E_UNASSIGNED = -1
integer(C_INT), public, bind(C) :: E_ALLOCATE = -2
integer(C_INT), public, bind(C) :: E_OUT_OF_BOUNDS = -3
integer(C_INT), public, bind(C) :: E_INVALID_SIZE = -4
integer(C_INT), public, bind(C) :: E_INVALID_ARGUMENT = -5
integer(C_INT), public, bind(C) :: E_INVALID_TYPE = -6
integer(C_INT), public, bind(C) :: E_INVALID_ID = -7
integer(C_INT), public, bind(C) :: E_GEOMETRY = -8
integer(C_INT), public, bind(C) :: E_DATA = -9
integer(C_INT), public, bind(C) :: E_PHYSICS = -10
integer(C_INT), public, bind(C, name='OPENMC_E_UNASSIGNED') :: E_UNASSIGNED = -1
integer(C_INT), public, bind(C, name='OPENMC_E_ALLOCATE') :: E_ALLOCATE = -2
integer(C_INT), public, bind(C, name='OPENMC_E_OUT_OF_BOUNDS') :: E_OUT_OF_BOUNDS = -3
integer(C_INT), public, bind(C, name='OPENMC_E_INVALID_SIZE') :: E_INVALID_SIZE = -4
integer(C_INT), public, bind(C, name='OPENMC_E_INVALID_ARGUMENT') :: E_INVALID_ARGUMENT = -5
integer(C_INT), public, bind(C, name='OPENMC_E_INVALID_TYPE') :: E_INVALID_TYPE = -6
integer(C_INT), public, bind(C, name='OPENMC_E_INVALID_ID') :: E_INVALID_ID = -7
integer(C_INT), public, bind(C, name='OPENMC_E_GEOMETRY') :: E_GEOMETRY = -8
integer(C_INT), public, bind(C, name='OPENMC_E_DATA') :: E_DATA = -9
integer(C_INT), public, bind(C, name='OPENMC_E_PHYSICS') :: E_PHYSICS = -10
! Warning codes
integer(C_INT), public, bind(C) :: E_WARNING = 1
integer(C_INT), public, bind(C, name='OPENMC_E_WARNING') :: E_WARNING = 1
! Error message
character(kind=C_CHAR), public, bind(C) :: openmc_err_msg(256)
@ -111,6 +111,14 @@ contains
end subroutine warning
subroutine warning_from_c(message, message_len) bind(C)
integer(C_INT), intent(in), value :: message_len
character(kind=C_CHAR), intent(in) :: message(message_len)
character(message_len+1) :: message_out
write(message_out, *) message
call warning(message_out)
end subroutine
!===============================================================================
! FATAL_ERROR alerts the user that an error has been encountered and displays a
! message about the particular problem. Errors are considered 'fatal' and hence

23
src/error.h
View file

@ -9,25 +9,38 @@
namespace openmc {
extern "C" void fatal_error_from_c(const char *message, int message_len);
extern "C" void fatal_error_from_c(const char* message, int message_len);
extern "C" void warning_from_c(const char* message, int message_len);
inline
void fatal_error(const char *message)
{
fatal_error_from_c(message, strlen(message));
}
inline
void fatal_error(const std::string &message)
{
fatal_error_from_c(message.c_str(), message.length());
}
inline
void fatal_error(const std::stringstream &message)
{
std::string out {message.str()};
fatal_error_from_c(out.c_str(), out.length());
fatal_error(message.str());
}
inline
void warning(const std::string& message)
{
warning_from_c(message.c_str(), message.length());
}
inline
void warning(const std::stringstream& message)
{
warning(message.str());
}
} // namespace openmc

10
src/finalize.cpp Normal file
View file

@ -0,0 +1,10 @@
#include "finalize.h"
#include "message_passing.h"
void openmc_free_bank()
{
#ifdef OPENMC_MPI
MPI_Type_free(&openmc::mpi::bank);
#endif
}

6
src/finalize.h Normal file
View file

@ -0,0 +1,6 @@
#ifndef FINALIZE_H
#define FINALIZE_H
extern "C" void openmc_free_bank();
#endif // FINALIZE_H

3
src/geometry.F90
View file

@ -355,9 +355,10 @@ contains
else
! Particle is outside the lattice.
if (lat % outer == NO_OUTER_UNIVERSE) then
call p % mark_as_lost("Particle " // trim(to_str(p %id)) &
call warning("Particle " // trim(to_str(p %id)) &
// " is outside lattice " // trim(to_str(lat % id)) &
// " but the lattice has no defined outer universe.")
found = .false.
return
else
p % coord(j + 1) % universe = lat % outer

14
src/geometry_header.F90
View file

@ -544,16 +544,12 @@ contains
c % type = FILL_MATERIAL
do i = 1, n
j = indices(i)
if (j == 0) then
c % material(i) = MATERIAL_VOID
if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then
c % material(i) = j
else
if (j >= 1 .and. j <= n_materials) then
c % material(i) = j
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index " // trim(to_str(j)) // " in the &
&materials array is out of bounds.")
end if
err = E_OUT_OF_BOUNDS
call set_errmsg("Index " // trim(to_str(j)) // " in the &
&materials array is out of bounds.")
end if
end do
case (FILL_UNIVERSE)

2803
src/hdf5_interface.F90
View file

File diff suppressed because it is too large Load diff

658
src/hdf5_interface.cpp Normal file
View file

@ -0,0 +1,658 @@
#include "hdf5_interface.h"
#include <array>
#include <cstring>
#include <sstream>
#include <string>
#include "hdf5.h"
#include "hdf5_hl.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#include "message_passing.h"
#endif
#include "error.h"
namespace openmc {
bool
attribute_exists(hid_t obj_id, const char* name)
{
htri_t out = H5Aexists_by_name(obj_id, ".", name, H5P_DEFAULT);
return out > 0;
}
size_t
attribute_typesize(hid_t obj_id, const char* name)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
hid_t filetype = H5Aget_type(attr);
size_t n = H5Tget_size(filetype);
H5Tclose(filetype);
H5Aclose(attr);
return n;
}
void
get_shape(hid_t obj_id, hsize_t* dims)
{
auto type = H5Iget_type(obj_id);
hid_t dspace;
if (type == H5I_DATASET) {
dspace = H5Dget_space(obj_id);
} else if (type == H5I_ATTR) {
dspace = H5Aget_space(obj_id);
}
H5Sget_simple_extent_dims(dspace, dims, nullptr);
H5Sclose(dspace);
}
void
get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
hid_t dspace = H5Aget_space(attr);
H5Sget_simple_extent_dims(dspace, dims, nullptr);
H5Sclose(dspace);
H5Aclose(attr);
}
hid_t
create_group(hid_t parent_id, char const *name)
{
hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (out < 0) {
std::stringstream err_msg;
err_msg << "Failed to create HDF5 group \"" << name << "\"";
fatal_error(err_msg);
}
return out;
}
hid_t
create_group(hid_t parent_id, const std::string &name)
{
return create_group(parent_id, name.c_str());
}
void
close_dataset(hid_t dataset_id)
{
if (H5Dclose(dataset_id) < 0) fatal_error("Failed to close dataset");
}
void
close_group(hid_t group_id)
{
if (H5Gclose(group_id) < 0) fatal_error("Failed to close group");
}
int
dataset_ndims(hid_t dset)
{
hid_t dspace = H5Dget_space(dset);
int ndims = H5Sget_simple_extent_ndims(dspace);
H5Sclose(dspace);
return ndims;
}
size_t
dataset_typesize(hid_t dset)
{
hid_t filetype = H5Dget_type(dset);
size_t n = H5Tget_size(filetype);
H5Tclose(filetype);
return n;
}
hid_t
file_open(const char* filename, char mode, bool parallel)
{
bool create;
unsigned int flags;
switch (mode) {
case 'r':
case 'a':
create = false;
flags = (mode == 'r' ? H5F_ACC_RDONLY : H5F_ACC_RDWR);
break;
case 'w':
case 'x':
create = true;
flags = (mode == 'x' ? H5F_ACC_EXCL : H5F_ACC_TRUNC);
break;
default:
std::stringstream err_msg;
err_msg << "Invalid file mode: " << mode;
fatal_error(err_msg);
}
hid_t plist = H5P_DEFAULT;
#ifdef PHDF5
if (parallel) {
// Setup file access property list with parallel I/O access
plist = H5Pcreate(H5P_FILE_ACCESS);
H5Pset_fapl_mpio(plist, openmc::mpi::intracomm, MPI_INFO_NULL);
}
#endif
// Open the file collectively
hid_t file_id;
if (create) {
file_id = H5Fcreate(filename, flags, H5P_DEFAULT, plist);
} else {
file_id = H5Fopen(filename, flags, plist);
}
if (file_id < 0) {
std::stringstream msg;
msg << "Failed to open HDF5 file with mode '" << mode << "': " << filename;
fatal_error(msg);
}
#ifdef PHDF5
// Close the property list
if (parallel) H5Pclose(plist);
#endif
return file_id;
}
hid_t
file_open(const std::string& filename, char mode, bool parallel=false)
{
file_open(filename.c_str(), mode, parallel);
}
void file_close(hid_t file_id)
{
H5Fclose(file_id);
}
void
get_name(hid_t obj_id, char* name)
{
size_t size = 1 + H5Iget_name(obj_id, nullptr, 0);
H5Iget_name(obj_id, name, size);
}
int get_num_datasets(hid_t group_id)
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get number of groups
H5O_info_t oinfo;
int ndatasets = 0;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type == H5O_TYPE_DATASET) ndatasets += 1;
}
return ndatasets;
}
int get_num_groups(hid_t group_id)
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get number of groups
H5O_info_t oinfo;
int ngroups = 0;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type == H5O_TYPE_GROUP) ngroups += 1;
}
return ngroups;
}
void
get_datasets(hid_t group_id, char* name[])
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get names
H5O_info_t oinfo;
hsize_t count = 0;
size_t size;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type != H5O_TYPE_DATASET) continue;
// Get size of name
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
i, nullptr, 0, H5P_DEFAULT);
// Read name
H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i,
name[count], size, H5P_DEFAULT);
count += 1;
}
}
void
get_groups(hid_t group_id, char* name[])
{
// Determine number of links in the group
H5G_info_t info;
H5Gget_info(group_id, &info);
// Iterate over links to get names
H5O_info_t oinfo;
hsize_t count = 0;
size_t size;
for (hsize_t i = 0; i < info.nlinks; ++i) {
// Determine type of object (and skip non-group)
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
H5P_DEFAULT);
if (oinfo.type != H5O_TYPE_GROUP) continue;
// Get size of name
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
i, nullptr, 0, H5P_DEFAULT);
// Read name
H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i,
name[count], size, H5P_DEFAULT);
count += 1;
}
}
bool
object_exists(hid_t object_id, const char* name)
{
htri_t out = H5LTpath_valid(object_id, name, true);
if (out < 0) {
std::stringstream err_msg;
err_msg << "Failed to check if object \"" << name << "\" exists.";
fatal_error(err_msg);
}
return (out > 0);
}
hid_t
open_dataset(hid_t group_id, const char* name)
{
if (object_exists(group_id, name)) {
return H5Dopen(group_id, name, H5P_DEFAULT);
} else {
std::stringstream err_msg;
err_msg << "Group \"" << name << "\" does not exist";
fatal_error(err_msg);
}
}
hid_t
open_group(hid_t group_id, const char* name)
{
if (object_exists(group_id, name)) {
return H5Gopen(group_id, name, H5P_DEFAULT);
} else {
std::stringstream err_msg;
err_msg << "Group \"" << name << "\" does not exist";
fatal_error(err_msg);
}
}
void
read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer)
{
hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT);
H5Aread(attr, mem_type_id, buffer);
H5Aclose(attr);
}
void
read_attr_double(hid_t obj_id, const char* name, double* buffer)
{
read_attr(obj_id, name, H5T_NATIVE_DOUBLE, buffer);
}
void
read_attr_int(hid_t obj_id, const char* name, int* buffer)
{
read_attr(obj_id, name, H5T_NATIVE_INT, buffer);
}
void
read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer)
{
// Create datatype for a string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen + 1);
// Read data into buffer
read_attr(obj_id, name, datatype, buffer);
// Free resources
H5Tclose(datatype);
}
void
read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
void* buffer, bool indep)
{
hid_t dset = obj_id;
if (name) dset = open_dataset(obj_id, name);
if (using_mpio_device(dset)) {
#ifdef PHDF5
// Set up collective vs independent I/O
auto data_xfer_mode = indep ? H5FD_MPIO_INDEPENDENT : H5FD_MPIO_COLLECTIVE;
// Create dataset transfer property list
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, data_xfer_mode);
// Read data
H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
H5Pclose(plist);
#endif
} else {
H5Dread(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
}
if (name) H5Dclose(dset);
}
void
read_double(hid_t obj_id, const char* name, double* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_DOUBLE, buffer, indep);
}
void
read_int(hid_t obj_id, const char* name, int* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_INT, buffer, indep);
}
void
read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep)
{
read_dataset(obj_id, name, H5T_NATIVE_LLONG, buffer, indep);
}
void
read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep)
{
// Create datatype for a string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen + 1);
// Read data into buffer
read_dataset(obj_id, name, datatype, buffer, indep);
// Free resources
H5Tclose(datatype);
}
void
read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep)
{
// Create compound datatype for complex numbers
struct complex_t {
double re;
double im;
};
complex_t tmp;
hid_t complex_id = H5Tcreate(H5T_COMPOUND, sizeof tmp);
H5Tinsert(complex_id, "r", HOFFSET(complex_t, re), H5T_NATIVE_DOUBLE);
H5Tinsert(complex_id, "i", HOFFSET(complex_t, im), H5T_NATIVE_DOUBLE);
// Read data
read_dataset(obj_id, name, complex_id, buffer, indep);
// Free resources
H5Tclose(complex_id);
}
void
read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
{
// Create dataspace for hyperslab in memory
hsize_t dims[] {n_filter, n_score, 3};
hsize_t start[] {0, 0, 1};
hsize_t count[] {n_filter, n_score, 2};
hid_t memspace = H5Screate_simple(3, dims, nullptr);
H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Create and write dataset
hid_t dset = H5Dopen(group_id, "results", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
// Free resources
H5Dclose(dset);
H5Sclose(memspace);
}
void
write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer)
{
// If array is given, create a simple dataspace. Otherwise, create a scalar
// datascape.
hid_t dspace;
if (ndim > 0) {
dspace = H5Screate_simple(ndim, dims, nullptr);
} else {
dspace = H5Screate(H5S_SCALAR);
}
// Create attribute and Write data
hid_t attr = H5Acreate(obj_id, name, mem_type_id, dspace,
H5P_DEFAULT, H5P_DEFAULT);
H5Awrite(attr, mem_type_id, buffer);
// Free resources
H5Aclose(attr);
H5Sclose(dspace);
}
void
write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
const double* buffer)
{
write_attr(obj_id, ndim, dims, name, H5T_NATIVE_DOUBLE, buffer);
}
void
write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
const int* buffer)
{
write_attr(obj_id, ndim, dims, name, H5T_NATIVE_INT, buffer);
}
void
write_attr_string(hid_t obj_id, const char* name, const char* buffer)
{
size_t n = strlen(buffer);
if (n > 0) {
// Set up appropriate datatype for a fixed-length string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, n);
write_attr(obj_id, 0, nullptr, name, datatype, buffer);
// Free resources
H5Tclose(datatype);
}
}
void
write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer, bool indep)
{
// If array is given, create a simple dataspace. Otherwise, create a scalar
// datascape.
hid_t dspace;
if (ndim > 0) {
dspace = H5Screate_simple(ndim, dims, nullptr);
} else {
dspace = H5Screate(H5S_SCALAR);
}
hid_t dset = H5Dcreate(group_id, name, mem_type_id, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (using_mpio_device(group_id)) {
#ifdef PHDF5
// Set up collective vs independent I/O
auto data_xfer_mode = indep ? H5FD_MPIO_INDEPENDENT : H5FD_MPIO_COLLECTIVE;
// Create dataset transfer property list
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, data_xfer_mode);
// Write data
H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, plist, buffer);
H5Pclose(plist);
#endif
} else {
H5Dwrite(dset, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
}
// Free resources
H5Dclose(dset);
H5Sclose(dspace);
}
void
write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const double* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_DOUBLE, buffer, indep);
}
void
write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const int* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_INT, buffer, indep);
}
void
write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
const long long* buffer, bool indep)
{
write_dataset(group_id, ndim, dims, name, H5T_NATIVE_LLONG, buffer, indep);
}
void
write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
const char* name, const char* buffer, bool indep)
{
if (slen > 0) {
// Set up appropriate datatype for a fixed-length string
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, slen);
write_dataset(group_id, ndim, dims, name, datatype, buffer, indep);
// Free resources
H5Tclose(datatype);
}
}
void
write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep)
{
write_string(group_id, 0, nullptr, buffer.length(), name, buffer.c_str(), indep);
}
void
write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results)
{
// Set dimensions of sum/sum_sq hyperslab to store
hsize_t count[] {n_filter, n_score, 2};
hid_t dspace = H5Screate_simple(3, count, nullptr);
// Set dimensions of results array
hsize_t dims[] {n_filter, n_score, 3};
hsize_t start[] {0, 0, 1};
hid_t memspace = H5Screate_simple(3, dims, nullptr);
H5Sselect_hyperslab(memspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Create and write dataset
hid_t dset = H5Dcreate(group_id, "results", H5T_NATIVE_DOUBLE, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dset, H5T_NATIVE_DOUBLE, memspace, H5S_ALL, H5P_DEFAULT, results);
// Free resources
H5Dclose(dset);
H5Sclose(memspace);
H5Sclose(dspace);
}
bool
using_mpio_device(hid_t obj_id)
{
// Determine file that this object is part of
hid_t file_id = H5Iget_file_id(obj_id);
// Get file access property list
hid_t fapl_id = H5Fget_access_plist(file_id);
// Get low-level driver identifier
hid_t driver = H5Pget_driver(fapl_id);
// Free resources
H5Pclose(fapl_id);
H5Fclose(file_id);
return driver == H5FD_MPIO;
}
} // namespace openmc

142
src/hdf5_interface.h
View file

@ -1,90 +1,96 @@
#ifndef HDF5_INTERFACE_H
#define HDF5_INTERFACE_H
#include <array>
#include <string.h>
#include "hdf5.h"
#include "hdf5_hl.h"
#include "error.h"
#include <array>
#include <string>
#include <sstream>
#include <complex.h>
namespace openmc {
extern "C" bool attribute_exists(hid_t obj_id, const char* name);
extern "C" size_t attribute_typesize(hid_t obj_id, const char* name);
extern "C" hid_t create_group(hid_t parent_id, const char* name);
hid_t create_group(hid_t parent_id, const std::string& name);
extern "C" void close_dataset(hid_t dataset_id);
extern "C" void close_group(hid_t group_id);
extern "C" int dataset_ndims(hid_t dset);
extern "C" size_t dataset_typesize(hid_t dset);
extern "C" hid_t file_open(const char* filename, char mode, bool parallel);
hid_t file_open(const std::string& filename, char mode, bool parallel);
extern "C" void file_close(hid_t file_id);
extern "C" void get_name(hid_t obj_id, char* name);
extern "C" int get_num_datasets(hid_t group_id);
extern "C" int get_num_groups(hid_t group_id);
extern "C" void get_datasets(hid_t group_id, char* name[]);
extern "C" void get_groups(hid_t group_id, char* name[]);
extern "C" void get_shape(hid_t obj_id, hsize_t* dims);
extern "C" void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims);
extern "C" bool object_exists(hid_t object_id, const char* name);
extern "C" hid_t open_dataset(hid_t group_id, const char* name);
extern "C" hid_t open_group(hid_t group_id, const char* name);
bool using_mpio_device(hid_t obj_id);
hid_t
create_group(hid_t parent_id, char const *name)
{
hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (out < 0) {
std::string err_msg{"Failed to create HDF5 group \""};
err_msg += name;
err_msg += "\"";
fatal_error(err_msg);
}
return out;
}
void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id,
const void* buffer);
extern "C" void read_attr_double(hid_t obj_id, const char* name, double* buffer);
extern "C" void read_attr_int(hid_t obj_id, const char* name, int* buffer);
extern "C" void read_attr_string(hid_t obj_id, const char* name, size_t slen,
char* buffer);
void read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
void* buffer, bool indep);
extern "C" void read_double(hid_t obj_id, const char* name, double* buffer,
bool indep);
extern "C" void read_int(hid_t obj_id, const char* name, int* buffer,
bool indep);
extern "C" void read_llong(hid_t obj_id, const char* name, long long* buffer,
bool indep);
extern "C" void read_string(hid_t obj_id, const char* name, size_t slen,
char* buffer, bool indep);
extern "C" void read_complex(hid_t obj_id, const char* name,
double _Complex* buffer, bool indep);
extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter,
hsize_t n_score, double* results);
void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer);
extern "C" void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
const char* name, const double* buffer);
extern "C" void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
const char* name, const int* buffer);
extern "C" void write_attr_string(hid_t obj_id, const char* name, const char* buffer);
hid_t
create_group(hid_t parent_id, const std::string &name)
{
return create_group(parent_id, name.c_str());
}
void write_dataset(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
hid_t mem_type_id, const void* buffer, bool indep);
extern "C" void write_double(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const double* buffer, bool indep);
extern "C" void write_int(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const int* buffer, bool indep);
extern "C" void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
const char* name, const long long* buffer, bool indep);
extern "C" void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
const char* name, char const* buffer, bool indep);
void write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep);
void
close_group(hid_t group_id)
{
herr_t err = H5Gclose(group_id);
if (err < 0) {
fatal_error("Failed to close HDF5 group");
}
}
extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
const double* results);
template<std::size_t array_len> void
write_double_1D(hid_t group_id, char const *name,
std::array<double, array_len> &buffer)
{
hsize_t dims[1]{array_len};
hid_t dataspace = H5Screate_simple(1, dims, NULL);
hid_t dataset = H5Dcreate(group_id, name, H5T_NATIVE_DOUBLE, dataspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&buffer[0]);
H5Sclose(dataspace);
H5Dclose(dataset);
}
void
write_string(hid_t group_id, char const *name, char const *buffer)
{
size_t buffer_len = strlen(buffer);
hid_t datatype = H5Tcopy(H5T_C_S1);
H5Tset_size(datatype, buffer_len);
hid_t dataspace = H5Screate(H5S_SCALAR);
hid_t dataset = H5Dcreate(group_id, name, datatype, dataspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, buffer);
H5Tclose(datatype);
H5Sclose(dataspace);
H5Dclose(dataset);
}
void
write_string(hid_t group_id, char const *name, const std::string &buffer)
{
write_string(group_id, name, buffer.c_str());
hsize_t dims[1] {array_len};
write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE,
buffer.data(), false);
}
} // namespace openmc

262
src/initialize.F90
View file

@ -1,37 +1,27 @@
module initialize
use, intrinsic :: ISO_C_BINDING, only: c_loc
use, intrinsic :: ISO_C_BINDING
use hdf5
#ifdef _OPENMP
use omp_lib
#endif
use bank_header, only: Bank
use constants
use set_header, only: SetInt
use error, only: fatal_error, warning, write_message
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
root_universe
use hdf5_interface, only: file_open, read_attribute, file_close, &
hdf5_bank_t, hdf5_integer8_t
use input_xml, only: read_input_xml
use material_header, only: Material
use message_passing
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: print_version, print_usage
use random_lcg, only: openmc_set_seed
use settings
#ifdef _OPENMP
use simulation_header, only: n_threads
#endif
use string, only: to_str, starts_with, ends_with, str_to_int
use tally_header, only: TallyObject
use tally_filter
use string, only: ends_with, to_f_string
use timer_header
implicit none
type(C_PTR), bind(C) :: openmc_path_input
type(C_PTR), bind(C) :: openmc_path_statepoint
type(C_PTR), bind(C) :: openmc_path_sourcepoint
type(C_PTR), bind(C) :: openmc_path_particle_restart
contains
!===============================================================================
@ -41,8 +31,9 @@ contains
! setting up timers, etc.
!===============================================================================
subroutine openmc_init(intracomm) bind(C)
function openmc_init_f(intracomm) result(err) bind(C)
integer, intent(in), optional :: intracomm ! MPI intracommunicator
integer(C_INT) :: err
#ifdef _OPENMP
character(MAX_WORD_LEN) :: envvar
@ -86,9 +77,6 @@ contains
end if
#endif
! Initialize HDF5 interface
call hdf5_initialize()
! Read command line arguments
call read_command_line()
@ -105,7 +93,8 @@ contains
! Stop initialization timer
call time_initialize%stop()
end subroutine openmc_init
err = 0
end function openmc_init_f
#ifdef OPENMC_MPI
!===============================================================================
@ -129,27 +118,11 @@ contains
integer :: bank_types(5) ! Datatypes
#endif
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
logical :: init_called
type(Bank) :: b
! Indicate that MPI is turned on
mpi_enabled = .true.
! Initialize MPI
call MPI_INITIALIZED(init_called, mpi_err)
if (.not. init_called) call MPI_INIT(mpi_err)
! Determine number of processors and rank of each processor
mpi_intracomm = intracomm
call MPI_COMM_SIZE(mpi_intracomm, n_procs, mpi_err)
call MPI_COMM_RANK(mpi_intracomm, rank, mpi_err)
! Determine master
if (rank == 0) then
master = .true.
else
master = .false.
end if
! ==========================================================================
! CREATE MPI_BANK TYPE
@ -174,203 +147,46 @@ contains
end subroutine initialize_mpi
#endif
!===============================================================================
! HDF5_INITIALIZE
!===============================================================================
subroutine hdf5_initialize()
type(Bank), target :: tmpb(2) ! temporary Bank
integer :: hdf5_err
integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals
integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates
! Initialize FORTRAN interface.
call h5open_f(hdf5_err)
! Create compound type for xyz and uvw
call h5tarray_create_f(H5T_NATIVE_DOUBLE, 1, dims, coordinates_t, hdf5_err)
! Create the compound datatype for Bank
call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(2))), hdf5_bank_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "wgt", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%wgt)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "xyz", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%xyz)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "uvw", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%uvw)), coordinates_t, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "E", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%E)), H5T_NATIVE_DOUBLE, hdf5_err)
call h5tinsert_f(hdf5_bank_t, "delayed_group", h5offsetof(c_loc(tmpb(1)), &
c_loc(tmpb(1)%delayed_group)), H5T_NATIVE_INTEGER, hdf5_err)
! Determine type for integer(8)
hdf5_integer8_t = h5kind_to_type(8, H5_INTEGER_KIND)
end subroutine hdf5_initialize
!===============================================================================
! READ_COMMAND_LINE reads all parameters from the command line
!===============================================================================
subroutine read_command_line()
! Arguments were already read on C++ side (initialize.cpp). Here we just
! convert the C-style strings to Fortran style
integer :: i ! loop index
integer :: argc ! number of command line arguments
integer :: last_flag ! index of last flag
character(MAX_WORD_LEN) :: filetype
integer(HID_T) :: file_id
character(MAX_WORD_LEN), allocatable :: argv(:) ! command line arguments
character(kind=C_CHAR), pointer :: string(:)
interface
function is_null(ptr) result(x) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: ptr
logical(C_BOOL) :: x
end function is_null
end interface
! Check number of command line arguments and allocate argv
argc = COMMAND_ARGUMENT_COUNT()
! Allocate and retrieve command arguments
allocate(argv(argc))
do i = 1, argc
call GET_COMMAND_ARGUMENT(i, argv(i))
end do
! Process command arguments
last_flag = 0
i = 1
do while (i <= argc)
! Check for flags
if (starts_with(argv(i), "-")) then
select case (argv(i))
case ('-p', '-plot', '--plot')
run_mode = MODE_PLOTTING
check_overlaps = .true.
case ('-n', '-particles', '--particles')
! Read number of particles per cycle
i = i + 1
n_particles = str_to_int(argv(i))
! Check that number specified was valid
if (n_particles == ERROR_INT) then
call fatal_error("Must specify integer after " // trim(argv(i-1)) &
&// " command-line flag.")
end if
case ('-r', '-restart', '--restart')
! Read path for state point/particle restart
i = i + 1
! Check what type of file this is
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
! Set path and flag for type of run
select case (trim(filetype))
case ('statepoint')
path_state_point = argv(i)
restart_run = .true.
case ('particle restart')
path_particle_restart = argv(i)
particle_restart_run = .true.
case default
call fatal_error("Unrecognized file after restart flag: " // filetype // ".")
end select
! If its a restart run check for additional source file
if (restart_run .and. i + 1 <= argc) then
! Increment arg
i = i + 1
! Check if it has extension we can read
if (ends_with(argv(i), '.h5')) then
! Check file type is a source file
file_id = file_open(argv(i), 'r', parallel=.true.)
call read_attribute(filetype, file_id, 'filetype')
call file_close(file_id)
if (filetype /= 'source') then
call fatal_error("Second file after restart flag must be a &
&source file")
end if
! It is a source file
path_source_point = argv(i)
else ! Different option is specified not a source file
! Source is in statepoint file
path_source_point = path_state_point
! Set argument back
i = i - 1
end if
else ! No command line arg after statepoint
! Source is assumed to be in statepoint file
path_source_point = path_state_point
end if
case ('-g', '-geometry-debug', '--geometry-debug')
check_overlaps = .true.
case ('-c', '--volume')
run_mode = MODE_VOLUME
case ('-s', '--threads')
! Read number of threads
i = i + 1
#ifdef _OPENMP
! Read and set number of OpenMP threads
n_threads = int(str_to_int(argv(i)), 4)
if (n_threads < 1) then
call fatal_error("Invalid number of threads specified on command &
&line.")
end if
call omp_set_num_threads(n_threads)
#else
if (master) call warning("Ignoring number of threads specified on &
&command line.")
#endif
case ('-?', '-h', '-help', '--help')
call print_usage()
stop
case ('-v', '-version', '--version')
call print_version()
stop
case ('-t', '-track', '--track')
write_all_tracks = .true.
case default
call fatal_error("Unknown command line option: " // argv(i))
end select
last_flag = i
end if
! Increment counter
i = i + 1
end do
! Determine directory where XML input files are
if (argc > 0 .and. last_flag < argc) then
path_input = argv(last_flag + 1)
if (.not. is_null(openmc_path_input)) then
call c_f_pointer(openmc_path_input, string, [255])
path_input = to_f_string(string)
else
path_input = ''
end if
! Add slash at end of directory if it isn't there
if (.not. ends_with(path_input, "/") .and. len_trim(path_input) > 0) then
path_input = trim(path_input) // "/"
if (.not. is_null(openmc_path_statepoint)) then
call c_f_pointer(openmc_path_statepoint, 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])
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])
path_particle_restart = to_f_string(string)
end if
! Free memory from argv
deallocate(argv)
! TODO: Check that directory exists
! Add slash at end of directory if it isn't there
if (len_trim(path_input) > 0 .and. .not. ends_with(path_input, "/")) then
path_input = trim(path_input) // "/"
end if
end subroutine read_command_line
end module initialize

229
src/initialize.cpp Normal file
View file

@ -0,0 +1,229 @@
#include "initialize.h"
#include <cstddef>
#include <cstring>
#include <iostream>
#include <sstream>
#include <string>
#include "error.h"
#include "hdf5_interface.h"
#include "message_passing.h"
#include "openmc.h"
#ifdef _OPENMP
#include "omp.h"
#endif
// data/functions from Fortran side
extern "C" bool openmc_check_overlaps;
extern "C" bool openmc_write_all_tracks;
extern "C" bool openmc_particle_restart_run;
extern "C" bool openmc_restart_run;
extern "C" void print_usage();
extern "C" void print_version();
// Paths to various files
extern "C" {
char* openmc_path_input;
char* openmc_path_statepoint;
char* openmc_path_sourcepoint;
char* openmc_path_particle_restart;
bool is_null(void* ptr) {return !ptr;}
}
int openmc_init(int argc, char* argv[], const void* intracomm)
{
#ifdef OPENMC_MPI
// Check if intracomm was passed
MPI_Comm comm;
if (intracomm) {
comm = *static_cast<const MPI_Comm *>(intracomm);
} else {
comm = MPI_COMM_WORLD;
}
// Initialize MPI for C++
openmc::initialize_mpi(comm);
#endif
// Parse command-line arguments
int err = openmc::parse_command_line(argc, argv);
if (err) return err;
// Continue with rest of initialization
#ifdef OPENMC_MPI
MPI_Fint fcomm = MPI_Comm_c2f(comm);
openmc_init_f(&fcomm);
#else
openmc_init_f(nullptr);
#endif
return 0;
}
namespace openmc {
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm)
{
openmc::mpi::intracomm = intracomm;
// Initialize MPI
int flag;
MPI_Initialized(&flag);
if (!flag) MPI_Init(nullptr, nullptr);
// Determine number of processes and rank for each
MPI_Comm_size(intracomm, &openmc::mpi::n_procs);
MPI_Comm_rank(intracomm, &openmc::mpi::rank);
// Set variable for Fortran side
openmc_n_procs = openmc::mpi::n_procs;
openmc_rank = openmc::mpi::rank;
openmc_master = (openmc::mpi::rank == 0);
// Create bank datatype
Bank b;
MPI_Aint disp[] {
offsetof(Bank, wgt),
offsetof(Bank, xyz),
offsetof(Bank, uvw),
offsetof(Bank, E),
offsetof(Bank, delayed_group)
};
int blocks[] {1, 3, 3, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
MPI_Type_create_struct(5, blocks, disp, types, &openmc::mpi::bank);
MPI_Type_commit(&openmc::mpi::bank);
}
#endif // OPENMC_MPI
inline bool ends_with(std::string const& value, std::string const& ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
int
parse_command_line(int argc, char* argv[])
{
char buffer[256]; // buffer for reading attribute
int last_flag = 0;
for (int i=1; i < argc; ++i) {
std::string arg {argv[i]};
if (arg[0] == '-') {
if (arg == "-p" || arg == "--plot") {
openmc_run_mode = RUN_MODE_PLOTTING;
openmc_check_overlaps = true;
} else if (arg == "-n" || arg == "--particles") {
i += 1;
n_particles = std::stoll(argv[i]);
} else if (arg == "-r" || arg == "--restart") {
i += 1;
// Check what type of file this is
hid_t file_id = file_open(argv[i], 'r', true);
size_t len = attribute_typesize(file_id, "filetype");
read_attr_string(file_id, "filetype", len, buffer);
file_close(file_id);
std::string filetype {buffer};
// Set path and flag for type of run
if (filetype == "statepoint") {
openmc_path_statepoint = argv[i];
openmc_restart_run = true;
} else if (filetype == "particle restart") {
openmc_path_particle_restart = argv[i];
openmc_particle_restart_run = true;
} else {
std::stringstream msg;
msg << "Unrecognized file after restart flag: " << filetype << ".";
strcpy(openmc_err_msg, msg.str().c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
// If its a restart run check for additional source file
if (openmc_restart_run && i + 1 < argc) {
// Check if it has extension we can read
if (ends_with(argv[i+1], ".h5")) {
// Check file type is a source file
file_id = file_open(argv[i+1], 'r', true);
len = attribute_typesize(file_id, "filetype");
read_attr_string(file_id, "filetype", len, buffer);
file_close(file_id);
if (filetype != "source") {
std::string msg {"Second file after restart flag must be a source file"};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
// It is a source file
openmc_path_sourcepoint = argv[i+1];
i += 1;
} else {
// Source is in statepoint file
openmc_path_sourcepoint = openmc_path_statepoint;
}
} else {
// Source is assumed to be in statepoint file
openmc_path_sourcepoint = openmc_path_statepoint;
}
} else if (arg == "-g" || arg == "--geometry-debug") {
openmc_check_overlaps = true;
} else if (arg == "-c" || arg == "--volume") {
openmc_run_mode = RUN_MODE_VOLUME;
} else if (arg == "-s" || arg == "--threads") {
// Read number of threads
i += 1;
#ifdef _OPENMP
// Read and set number of OpenMP threads
openmc_n_threads = std::stoi(argv[i]);
if (openmc_n_threads < 1) {
std::string msg {"Number of threads must be positive."};
strcpy(openmc_err_msg, msg.c_str());
return OPENMC_E_INVALID_ARGUMENT;
}
omp_set_num_threads(openmc_n_threads);
#else
if (openmc_master)
warning("Ignoring number of threads specified on command line.");
#endif
} else if (arg == "-?" || arg == "-h" || arg == "--help") {
print_usage();
return OPENMC_E_UNASSIGNED;
} else if (arg == "-v" || arg == "--version") {
print_version();
return OPENMC_E_UNASSIGNED;
} else if (arg == "-t" || arg == "--track") {
openmc_write_all_tracks = true;
} else {
std::cerr << "Unknown option: " << argv[i] << '\n';
print_usage();
return OPENMC_E_UNASSIGNED;
}
last_flag = i;
}
}
// Determine directory where XML input files are
if (argc > 1 && last_flag < argc) openmc_path_input = argv[last_flag + 1];
return 0;
}
} // namespace openmc

20
src/initialize.h Normal file
View file

@ -0,0 +1,20 @@
#ifndef INITIALIZE_H
#define INITIALIZE_H
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
extern "C" void print_usage();
extern "C" void print_version();
namespace openmc {
int parse_command_line(int argc, char* argv[]);
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm);
#endif
}
#endif // INITIALIZE_H

240
src/input_xml.F90
View file

@ -2194,26 +2194,20 @@ contains
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: k ! another loop index
integer :: l ! another loop index
integer :: filter_id ! user-specified identifier for filter
integer :: i_filt ! index in filters array
integer :: i_elem ! index of entry in dictionary
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
integer :: n_filter ! number of filters
integer :: n_new ! number of new scores to add based on Yn/Pn tally
integer :: n_scores ! number of tot scores after adjusting for Yn/Pn tally
integer :: n_bins ! total new bins for this score
integer :: n_scores ! number of scores
integer :: n_user_trig ! number of user-specified tally triggers
integer :: trig_ind ! index of triggers array for each tally
integer :: user_trig_ind ! index of user-specified triggers for each tally
integer :: i_start, i_end
integer(C_INT) :: err
real(8) :: threshold ! trigger convergence threshold
integer :: n_order ! moment order requested
integer :: n_order_pos ! oosition of Scattering order in score name string
integer :: MT ! user-specified MT for score
integer :: imomstr ! Index of MOMENT_STRS & MOMENT_N_STRS
logical :: file_exists ! does tallies.xml file exist?
integer, allocatable :: temp_filter(:) ! temporary filter indices
character(MAX_LINE_LEN) :: filename
@ -2541,107 +2535,22 @@ contains
allocate(sarray(n_words))
call get_node_array(node_tal, "scores", sarray)
! Before we can allocate storage for scores, we must determine the
! number of additional scores required due to the moment scores
! (i.e., scatter-p#, flux-y#)
n_new = 0
! Append the score to the list of possible trigger scores
do j = 1, n_words
sarray(j) = to_lower(sarray(j))
! Find if scores(j) is of the form 'moment-p' or 'moment-y' present in
! MOMENT_STRS(:)
! If so, check the order, store if OK, then reset the number to 'n'
score_name = trim(sarray(j))
! Append the score to the list of possible trigger scores
if (trigger_on) call trigger_scores % set(trim(score_name), j)
do imomstr = 1, size(MOMENT_STRS)
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
! The above scheme will essentially take the absolute value
if (master) call warning("Invalid scattering order of " &
// trim(to_str(n_order)) // " requested. Setting to the &
&maximum permissible value, " &
// trim(to_str(MAX_ANG_ORDER)))
n_order = MAX_ANG_ORDER
sarray(j) = trim(MOMENT_STRS(imomstr)) &
// trim(to_str(MAX_ANG_ORDER))
end if
! Find total number of bins for this case
if (imomstr >= YN_LOC) then
n_bins = (n_order + 1)**2
else
n_bins = n_order + 1
end if
! We subtract one since n_words already included
n_new = n_new + n_bins - 1
exit
end if
end do
end do
n_scores = n_words + n_new
n_scores = n_words
! Allocate score storage accordingly
allocate(t % score_bins(n_scores))
allocate(t % moment_order(n_scores))
t % moment_order = 0
j = 0
do l = 1, n_words
j = j + 1
! Get the input string in scores(l) but if score is one of the moment
! scores then strip off the n and store it as an integer to be used
! later. Then perform the select case on this modified (number
! removed) string
n_order = -1
score_name = sarray(l)
do imomstr = 1, size(MOMENT_STRS)
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
! The above scheme will essentially take the absolute value
n_order = MAX_ANG_ORDER
end if
score_name = trim(MOMENT_STRS(imomstr)) // "n"
! Find total number of bins for this case
if (imomstr >= YN_LOC) then
n_bins = (n_order + 1)**2
else
n_bins = n_order + 1
end if
exit
end if
end do
! Now check the Moment_N_Strs, but only if we werent successful above
if (imomstr > size(MOMENT_STRS)) then
do imomstr = 1, size(MOMENT_N_STRS)
if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
! The above scheme will essentially take the absolute value
if (master) call warning("Invalid scattering order of " &
// trim(to_str(n_order)) // " requested. Setting to &
&the maximum permissible value, " &
// trim(to_str(MAX_ANG_ORDER)))
n_order = MAX_ANG_ORDER
end if
score_name = trim(MOMENT_N_STRS(imomstr)) // "n"
exit
end if
end do
end if
! Check the validity of the scores and their filters
do j = 1, n_scores
score_name = sarray(j)
! Check if delayed group filter is used with any score besides
! delayed-nu-fission or decay-rate
@ -2652,23 +2561,6 @@ contains
&delayedgroup filter.")
end if
! Check to see if the mu filter is applied and if that makes sense.
if ((.not. starts_with(score_name,'scatter')) .and. &
(.not. starts_with(score_name,'nu-scatter'))) then
if (t % find_filter(FILTER_MU) > 0) then
call fatal_error("Cannot tally " // trim(score_name) //" with a &
&change of angle (mu) filter.")
end if
! Also check to see if this is a legendre expansion or not.
! If so, we can accept this score and filter combo for p0, but not
! elsewhere.
else if (n_order > 0) then
if (t % find_filter(FILTER_MU) > 0) then
call fatal_error("Cannot tally " // trim(score_name) //" with a &
&change of angle (mu) filter unless order is 0.")
end if
end if
select case (trim(score_name))
case ('flux')
! Prohibit user from tallying flux for an individual nuclide
@ -2683,22 +2575,6 @@ contains
&filter.")
end if
case ('flux-yn')
! Prohibit user from tallying flux for an individual nuclide
if (.not. (t % n_nuclide_bins == 1 .and. &
t % nuclide_bins(1) == -1)) then
call fatal_error("Cannot tally flux for an individual nuclide.")
end if
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally flux with an outgoing energy &
&filter.")
end if
t % score_bins(j : j + n_bins - 1) = SCORE_FLUX_YN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case ('total', '(n,total)')
t % score_bins(j) = SCORE_TOTAL
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
@ -2706,18 +2582,13 @@ contains
&outgoing energy filter.")
end if
case ('total-yn')
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally total reaction rate with an &
&outgoing energy filter.")
end if
t % score_bins(j : j + n_bins - 1) = SCORE_TOTAL_YN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case ('scatter')
t % score_bins(j) = SCORE_SCATTER
if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
t % find_filter(FILTER_LEGENDRE) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('nu-scatter')
t % score_bins(j) = SCORE_NU_SCATTER
@ -2727,53 +2598,14 @@ contains
! necessary)
if (run_CE) then
t % estimator = ESTIMATOR_ANALOG
else
if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
t % find_filter(FILTER_LEGENDRE) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
end if
case ('scatter-n')
t % score_bins(j) = SCORE_SCATTER_N
t % moment_order(j) = n_order
t % estimator = ESTIMATOR_ANALOG
case ('nu-scatter-n')
t % score_bins(j) = SCORE_NU_SCATTER_N
t % moment_order(j) = n_order
t % estimator = ESTIMATOR_ANALOG
case ('scatter-pn')
t % estimator = ESTIMATOR_ANALOG
! Setup P0:Pn
t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_PN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case ('nu-scatter-pn')
t % estimator = ESTIMATOR_ANALOG
! Setup P0:Pn
t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_PN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case ('scatter-yn')
t % estimator = ESTIMATOR_ANALOG
! Setup P0:Pn
t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_YN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case ('nu-scatter-yn')
t % estimator = ESTIMATOR_ANALOG
! Setup P0:Pn
t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_YN
t % moment_order(j : j + n_bins - 1) = n_order
j = j + n_bins - 1
case('transport')
call fatal_error("Transport score no longer supported for tallies, &
&please remove")
case ('n1n')
call fatal_error("n1n score no longer supported for tallies, &
&please remove")
case ('n2n', '(n,2n)')
t % score_bins(j) = N_2N
t % depletion_rx = .true.
@ -2937,6 +2769,14 @@ contains
t % score_bins(j) = N_DA
case default
! First look for deprecated scores
if (starts_with(trim(score_name), 'scatter-') .or. &
starts_with(trim(score_name), 'nu-scatter-') .or. &
starts_with(trim(score_name), 'total-y') .or. &
starts_with(trim(score_name), 'flux-y')) then
call fatal_error(trim(score_name) // " is no longer available.")
end if
! Assume that user has specified an MT number
MT = int(str_to_int(score_name))
@ -2945,14 +2785,12 @@ contains
if (MT > 1) then
t % score_bins(j) = MT
else
call fatal_error("Invalid MT on <scores>: " &
// trim(sarray(l)))
call fatal_error("Invalid MT on <scores>: " // trim(score_name))
end if
else
! Specified score was not an integer
call fatal_error("Unknown scoring function: " &
// trim(sarray(l)))
call fatal_error("Unknown scoring function: " // trim(score_name))
end if
end select
@ -2967,35 +2805,19 @@ contains
end do
t % n_score_bins = n_scores
t % n_user_score_bins = n_words
! Deallocate temporary string array of scores
deallocate(sarray)
! Check that no duplicate scores exist
j = 1
do while (j < n_scores)
! Determine number of bins for scores with expansions
n_order = t % moment_order(j)
select case (t % score_bins(j))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
n_bins = n_order + 1
case (SCORE_FLUX_YN, SCORE_TOTAL_YN, SCORE_SCATTER_YN, &
SCORE_NU_SCATTER_YN)
n_bins = (n_order + 1)**2
case default
n_bins = 1
end select
do k = j + n_bins, n_scores
if (t % score_bins(j) == t % score_bins(k) .and. &
t % moment_order(j) == t % moment_order(k)) then
do j = 1, n_scores - 1
do k = j + 1, n_scores
if (t % score_bins(j) == t % score_bins(k)) then
call fatal_error("Duplicate score of type '" // trim(&
reaction_name(t % score_bins(j))) // "' found in tally " &
// trim(to_str(t % id)))
end if
end do
j = j + n_bins
end do
else
call fatal_error("No <scores> specified on tally " &

48
src/main.F90
View file

@ -1,48 +0,0 @@
program main
use constants
use message_passing
use openmc_api, only: openmc_init, openmc_finalize, openmc_run, &
openmc_plot_geometry, openmc_calculate_volumes
use particle_restart, only: run_particle_restart
use settings, only: run_mode
implicit none
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
! Initialize run -- when run with MPI, pass communicator
#ifdef OPENMC_MPI
#ifdef OPENMC_MPIF08
call openmc_init(MPI_COMM_WORLD % MPI_VAL)
#else
call openmc_init(MPI_COMM_WORLD)
#endif
#else
call openmc_init()
#endif
! start problem based on mode
select case (run_mode)
case (MODE_FIXEDSOURCE, MODE_EIGENVALUE)
call openmc_run()
case (MODE_PLOTTING)
call openmc_plot_geometry()
case (MODE_PARTICLE)
if (master) call run_particle_restart()
case (MODE_VOLUME)
call openmc_calculate_volumes()
end select
! finalize run
call openmc_finalize()
#ifdef OPENMC_MPI
! If MPI is in use and enabled, terminate it
call MPI_FINALIZE(mpi_err)
#endif
end program main

51
src/main.cpp Normal file
View file

@ -0,0 +1,51 @@
#include "error.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "openmc.h"
int main(int argc, char* argv[]) {
int err;
// Initialize run -- when run with MPI, pass communicator
#ifdef OPENMC_MPI
MPI_Comm world {MPI_COMM_WORLD};
err = openmc_init(argc, argv, &world);
#else
err = openmc_init(argc, argv, nullptr);
#endif
if (err == -1) {
// This happens for the -h and -v flags
return 0;
} else if (err) {
openmc::fatal_error(openmc_err_msg);
}
// start problem based on mode
switch (openmc_run_mode) {
case RUN_MODE_FIXEDSOURCE:
case RUN_MODE_EIGENVALUE:
err = openmc_run();
break;
case RUN_MODE_PLOTTING:
err = openmc_plot_geometry();
break;
case RUN_MODE_PARTICLE:
if (openmc_master) err = openmc_particle_restart();
break;
case RUN_MODE_VOLUME:
err = openmc_calculate_volumes();
break;
}
if (err) openmc::fatal_error(openmc_err_msg);
// Finalize and free up memory
err = openmc_finalize();
if (err) openmc::fatal_error(openmc_err_msg);
// If MPI is in use and enabled, terminate it
#ifdef OPENMC_MPI
MPI_Finalize();
#endif
}

13
src/math.F90
View file

@ -598,15 +598,8 @@ contains
real(8) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is
! easier to work with
real(8) :: k1, k2, k3, k4 ! Variables for R_m_n calculation
real(8) :: sqrt_norm ! normalization for radial moments
integer :: i,p,q ! Loop counters
real(8), parameter :: SQRT_N_1(0:10) = [&
sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), &
sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), &
sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)]
real(8), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8)
! n == radial degree
! m == azimuthal frequency
@ -670,11 +663,11 @@ contains
do p = 0, n
do q = -p, p, 2
if (q < 0) then
zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q)) * SQRT_2N_2(p)
zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q))
else if (q == 0) then
zn(i) = zn_mat(p+1, q+1) * SQRT_N_1(p)
zn(i) = zn_mat(p+1, q+1)
else
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) * SQRT_2N_2(p)
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1)
end if
i = i + 1
end do

199
src/mesh_header.F90
View file

@ -2,11 +2,9 @@ module mesh_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use dict_header, only: DictIntInt
use error, only: warning, fatal_error
use error
use hdf5_interface
use string, only: to_str, to_lower
use xml_interface
@ -15,6 +13,13 @@ module mesh_header
private
public :: free_memory_mesh
public :: openmc_extend_meshes
public :: openmc_get_mesh_index
public :: openmc_mesh_get_id
public :: openmc_mesh_get_dimension
public :: openmc_mesh_get_params
public :: openmc_mesh_set_id
public :: openmc_mesh_set_dimension
public :: openmc_mesh_set_params
!===============================================================================
! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
@ -23,13 +28,13 @@ module mesh_header
type, public :: RegularMesh
integer :: id = -1 ! user-specified id
integer :: type ! rectangular, hexagonal
integer :: n_dimension ! rank of mesh
integer :: type = MESH_REGULAR ! rectangular, hexagonal
integer(C_INT) :: n_dimension ! rank of mesh
real(8) :: volume_frac ! volume fraction of each cell
integer, allocatable :: dimension(:) ! number of cells in each direction
real(8), allocatable :: lower_left(:) ! lower-left corner of mesh
real(8), allocatable :: upper_right(:) ! upper-right corner of mesh
real(8), allocatable :: width(:) ! width of each mesh cell
integer(C_INT), allocatable :: dimension(:) ! number of cells in each direction
real(C_DOUBLE), allocatable :: lower_left(:) ! lower-left corner of mesh
real(C_DOUBLE), allocatable :: upper_right(:) ! upper-right corner of mesh
real(C_DOUBLE), allocatable :: width(:) ! width of each mesh cell
contains
procedure :: from_xml => regular_from_xml
procedure :: get_bin => regular_get_bin
@ -592,4 +597,180 @@ contains
err = 0
end function openmc_extend_meshes
function openmc_get_mesh_index(id, index) result(err) bind(C)
! Return the index in the meshes array of a mesh with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(meshes)) then
if (mesh_dict % has(id)) then
index = mesh_dict % get(id)
err = 0
else
err = E_INVALID_ID
call set_errmsg("No mesh exists with ID=" // trim(to_str(id)) // ".")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory has not been allocated for meshes.")
end if
end function openmc_get_mesh_index
function openmc_mesh_get_id(index, id) result(err) bind(C)
! Return the ID of a mesh
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(meshes)) then
id = meshes(index) % id
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_get_id
function openmc_mesh_set_id(index, id) result(err) bind(C)
! Set the ID of a mesh
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= n_meshes) then
meshes(index) % id = id
call mesh_dict % set(id, index)
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_set_id
function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C)
! Get the dimension of a mesh
integer(C_INT32_T), value, intent(in) :: index
type(C_PTR), intent(out) :: dims
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
if (index >= 1 .and. index <= n_meshes) then
dims = C_LOC(meshes(index) % dimension)
n = meshes(index) % n_dimension
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_get_dimension
function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C)
! Set the dimension of a mesh
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
integer(C_INT), intent(in) :: dims(n)
integer(C_INT) :: err
if (index >= 1 .and. index <= n_meshes) then
associate (m => meshes(index))
if (allocated(m % dimension)) deallocate (m % dimension)
if (allocated(m % lower_left)) deallocate (m % lower_left)
if (allocated(m % upper_right)) deallocate (m % upper_right)
if (allocated(m % width)) deallocate (m % width)
m % n_dimension = n
allocate(m % dimension(n))
allocate(m % lower_left(n))
allocate(m % upper_right(n))
allocate(m % width(n))
! Copy dimension
m % dimension(:) = dims
end associate
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_set_dimension
function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C)
! Get the mesh parameters
integer(C_INT32_T), value, intent(in) :: index
type(C_PTR), intent(out) :: ll
type(C_PTR), intent(out) :: ur
type(C_PTR), intent(out) :: width
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= n_meshes) then
associate (m => meshes(index))
if (allocated(m % lower_left)) then
ll = C_LOC(m % lower_left(1))
ur = C_LOC(m % upper_right(1))
width = C_LOC(m % width(1))
n = m % n_dimension
else
err = E_ALLOCATE
call set_errmsg("Mesh parameters have not been set.")
end if
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_get_params
function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C)
! Set the mesh parameters
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in), optional :: ll(n)
real(C_DOUBLE), intent(in), optional :: ur(n)
real(C_DOUBLE), intent(in), optional :: width(n)
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= n_meshes) then
associate (m => meshes(index))
if (allocated(m % lower_left)) deallocate (m % lower_left)
if (allocated(m % upper_right)) deallocate (m % upper_right)
if (allocated(m % width)) deallocate (m % width)
allocate(m % lower_left(n))
allocate(m % upper_right(n))
allocate(m % width(n))
if (present(ll) .and. present(ur)) then
m % lower_left(:) = ll
m % upper_right(:) = ur
m % width(:) = (ur - ll) / m % dimension
elseif (present(ll) .and. present(width)) then
m % lower_left(:) = ll
m % width(:) = width
m % upper_right(:) = ll + width * m % dimension
elseif (present(ur) .and. present(width)) then
m % upper_right(:) = ur
m % width(:) = width
m % lower_left(:) = ur - width * m % dimension
else
err = E_INVALID_ARGUMENT
call set_errmsg("At least two parameters must be specified.")
end if
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in meshes array is out of bounds.")
end if
end function openmc_mesh_set_params
end module mesh_header

8
src/message_passing.F90
View file

@ -1,5 +1,7 @@
module message_passing
use, intrinsic :: ISO_C_BINDING
#ifdef OPENMC_MPI
#ifdef OPENMC_MPIF08
use mpi_f08
@ -12,9 +14,9 @@ module message_passing
! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial
! run. In this case, these variables are still used at times.
integer :: n_procs = 1 ! number of processes
integer :: rank = 0 ! rank of process
logical :: master = .true. ! master process?
integer(C_INT), bind(C, name='openmc_n_procs') :: n_procs = 1 ! number of processes
integer(C_INT), bind(C, name='openmc_rank') :: rank = 0 ! rank of process
logical(C_BOOL), bind(C, name='openmc_master') :: master = .true. ! master process?
logical :: mpi_enabled = .false. ! is MPI in use and initialized?
#ifdef OPENMC_MPIF08
type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank

15
src/message_passing.cpp Normal file
View file

@ -0,0 +1,15 @@
#include "message_passing.h"
namespace openmc {
namespace mpi {
int rank {0};
int n_procs {1};
#ifdef OPENMC_MPI
MPI_Comm intracomm;
MPI_Datatype bank;
#endif
} // namespace mpi
} // namespace openmc

22
src/message_passing.h Normal file
View file

@ -0,0 +1,22 @@
#ifndef MESSAGE_PASSING_H
#define MESSAGE_PASSING_H
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
namespace openmc {
namespace mpi {
extern int rank;
extern int n_procs;
#ifdef OPENMC_MPI
extern MPI_Datatype bank;
extern MPI_Comm intracomm;
#endif
} // namespace mpi
} // namespace openmc
#endif // MESSAGE_PASSING_H

4
src/mgxs_data.F90
View file

@ -132,6 +132,8 @@ contains
! Add name to dictionary
call already_read % add(name)
call close_group(xsdata_group)
end if
end do NUCLIDE_LOOP
end do MATERIAL_LOOP
@ -159,6 +161,8 @@ contains
end do NUCLIDE_LOOP2
end do MATERIAL_LOOP3
call file_close(file_id)
end subroutine read_mgxs
!===============================================================================

36
src/mgxs_header.F90
View file

@ -3,10 +3,8 @@ module mgxs_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: find, sort
use constants, only: MAX_WORD_LEN, ZERO, ONE, TWO, PI
use constants, only: MAX_WORD_LEN, ZERO, ONE, TWO, PI, MACROSCOPIC_AWR
use error, only: fatal_error
use hdf5_interface
use material_header, only: material
@ -248,7 +246,6 @@ contains
type(VectorInt), intent(out) :: temps_to_read ! Temperatures to read
integer, intent(out) :: order_dim ! Scattering data order size
integer(SIZE_T) :: name_len
integer(HID_T) :: kT_group
character(MAX_WORD_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
@ -259,8 +256,7 @@ contains
integer :: ipol, iazi
! Get name of dataset from group
name_len = len(this % name)
this % name = get_name(xs_id, name_len)
this % name = get_name(xs_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -268,7 +264,7 @@ contains
if (attribute_exists(xs_id, "atomic_weight_ratio")) then
call read_attribute(this % awr, xs_id, "atomic_weight_ratio")
else
this % awr = -ONE
this % awr = MACROSCOPIC_AWR
end if
! Determine temperatures available
@ -396,9 +392,9 @@ contains
! Store the dimensionality of the data in order_dim.
! For Legendre data, we usually refer to it as Pn where n is the order.
! However Pn has n+1 sets of points (since you need to
! the count the P0 moment). Adjust for that. Histogram and Tabular
! formats dont need this adjustment.
! However Pn has n+1 sets of points (since you need to count the P0
! moment). Adjust for that. Histogram and Tabular formats dont need this
! adjustment.
if (this % scatter_format == ANGLE_LEGENDRE) then
order_dim = order_dim + 1
else
@ -551,6 +547,8 @@ contains
else
call fatal_error("beta must be provided as a 1D or 2D array")
end if
call close_dataset(xsdata)
else
temp_beta = ZERO
end if
@ -681,6 +679,8 @@ contains
call fatal_error("nu-fission must be provided as a 1D or 2D &
&array")
end if
call close_dataset(xsdata)
end if
! If chi-prompt provided, set chi-prompt
@ -784,6 +784,8 @@ contains
call fatal_error("chi-delayed must be provided as a 1D or 2D &
&array")
end if
call close_dataset(xsdata)
end if
! If prompt-nu-fission present, set prompt-nu-fission
@ -840,6 +842,8 @@ contains
call fatal_error("prompt-nu-fission must be provided as a 1D &
&or 2D array")
end if
call close_dataset(xsdata)
end if
! If delayed-nu-fission provided, set delayed-nu-fission. If
@ -965,6 +969,8 @@ contains
call fatal_error("delayed-nu-fission must be provided as a &
&1D, 2D, or 3D array")
end if
call close_dataset(xsdata)
end if
! Deallocate temporary beta array
@ -1348,6 +1354,8 @@ contains
else
call fatal_error("beta must be provided as a 3D or 4D array")
end if
call close_dataset(xsdata)
else
temp_beta = ZERO
end if
@ -1519,6 +1527,8 @@ contains
call fatal_error("nu-fission must be provided as a 3D or &
&4D array")
end if
call close_dataset(xsdata)
end if
! If chi-prompt provided, set chi-prompt
@ -1653,6 +1663,8 @@ contains
call fatal_error("chi-delayed must be provided as a 3D or 4D &
&array")
end if
call close_dataset(xsdata)
end if
! If prompt-nu-fission present, set prompt-nu-fission
@ -1722,6 +1734,8 @@ contains
call fatal_error("prompt-nu-fission must be provided as a 3D &
&or 4D array")
end if
call close_dataset(xsdata)
end if
! If delayed-nu-fission provided, set delayed-nu-fission. If
@ -1870,6 +1884,8 @@ contains
call fatal_error("delayed-nu-fission must be provided as a &
&3D, 4D, or 5D array")
end if
call close_dataset(xsdata)
end if
! Deallocate temporary beta array

2
src/multipole_header.F90
View file

@ -1,7 +1,5 @@
module multipole_header
use hdf5
use constants
use dict_header, only: DictIntInt
use error, only: fatal_error

8
src/nuclide_header.F90
View file

@ -3,8 +3,6 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: sort, find, binary_search
use constants
use dict_header, only: DictIntInt, DictCharInt
@ -277,7 +275,7 @@ contains
integer, intent(inout) :: method
real(8), intent(in) :: tolerance
real(8), intent(in) :: minmax(2) ! range of temperatures
logical, intent(in) :: master ! if this is the master proc
logical(C_BOOL), intent(in) :: master ! if this is the master proc
integer, intent(in) :: i_nuclide ! Nuclide index in nuclides
integer :: i
@ -292,7 +290,6 @@ contains
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp_str
@ -306,8 +303,7 @@ contains
type(VectorInt) :: index_inelastic_scatter
! Get name of nuclide from group
name_len = len(this % name)
this % name = get_name(group_id, name_len)
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))

75
src/output.F90
View file

@ -162,7 +162,7 @@ contains
! information
!===============================================================================
subroutine print_version()
subroutine print_version() bind(C)
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I2,".",I1)') &
@ -182,7 +182,7 @@ contains
! PRINT_USAGE displays information about command line usage of OpenMC
!===============================================================================
subroutine print_usage()
subroutine print_usage() bind(C)
if (master) then
write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]'
@ -665,14 +665,11 @@ contains
integer :: j ! level in tally hierarchy
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: l ! loop index for user scores
integer :: h ! loop index for tally filters
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
integer :: i_nuclide ! index in nuclides array
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer :: unit_tally ! tallies.out file unit
integer :: nr ! number of realizations
real(8) :: t_value ! t-values for confidence intervals
@ -699,14 +696,6 @@ contains
score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate"
score_names(abs(SCORE_EVENTS)) = "Events"
score_names(abs(SCORE_FLUX_YN)) = "Flux Moment"
score_names(abs(SCORE_TOTAL_YN)) = "Total Reaction Rate Moment"
score_names(abs(SCORE_SCATTER_N)) = "Scattering Rate Moment"
score_names(abs(SCORE_SCATTER_PN)) = "Scattering Rate Moment"
score_names(abs(SCORE_SCATTER_YN)) = "Scattering Rate Moment"
score_names(abs(SCORE_NU_SCATTER_N)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment"
score_names(abs(SCORE_DECAY_RATE)) = "Decay Rate"
score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
score_names(abs(SCORE_PROMPT_NU_FISSION)) = "Prompt-Nu-Fission Rate"
@ -860,60 +849,20 @@ contains
end if
indent = indent + 2
k = 0
do l = 1, t % n_user_score_bins
k = k + 1
do k = 1, t % n_score_bins
score_index = score_index + 1
associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
select case(t % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // &
score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, to_str(x(1)), &
trim(to_str(t_value * x(2)))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
score_index = score_index + 1
score_name = 'P' // trim(to_str(n_order)) // " " //&
score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(x(1)), trim(to_str(t_value * x(2)))
end do
k = k + t % moment_order(k)
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
do nm_order = -n_order, n_order
score_index = score_index + 1
score_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order)) // " " &
// score_names(abs(t % score_bins(k)))
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(x(1)), trim(to_str(t_value * x(2)))
end do
end do
k = k + (t % moment_order(k) + 1)**2 - 1
case default
if (t % score_bins(k) > 0) then
score_name = reaction_name(t % score_bins(k))
else
score_name = score_names(abs(t % score_bins(k)))
end if
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(x(1)), trim(to_str(t_value * x(2)))
end select
if (t % score_bins(k) > 0) then
score_name = reaction_name(t % score_bins(k))
else
score_name = score_names(abs(t % score_bins(k)))
end if
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(x(1)), trim(to_str(t_value * x(2)))
end associate
end do
indent = indent - 2

4
src/particle_header.F90
View file

@ -1,7 +1,5 @@
module particle_header
use hdf5, only: HID_T
use bank_header, only: Bank, source_bank
use constants
use error, only: fatal_error, warning
@ -312,7 +310,7 @@ contains
!$omp critical (WriteParticleRestart)
! Create file
file_id = file_create(filename)
file_id = file_open(filename, 'w')
associate (src => source_bank(current_work))
! Write filetype and version info

15
src/particle_restart.F90
View file

@ -5,7 +5,7 @@ module particle_restart
use bank_header, only: Bank
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_header, only: energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use output, only: print_particle
@ -16,24 +16,25 @@ module particle_restart
use tally_header, only: n_tallies
use tracking, only: transport
use hdf5, only: HID_T
implicit none
private
public :: run_particle_restart
public :: openmc_particle_restart
contains
!===============================================================================
! RUN_PARTICLE_RESTART is the main routine that runs the particle restart
! OPENMC_PARTICLE_RESTART is the main routine that runs the particle restart
!===============================================================================
subroutine run_particle_restart()
function openmc_particle_restart() result(err) bind(C)
integer(C_INT) :: err
integer(8) :: particle_seed
integer :: previous_run_mode
type(Particle) :: p
err = 0
! Set verbosity high
verbosity = 10
@ -66,7 +67,7 @@ contains
deallocate(micro_xs)
end subroutine run_particle_restart
end function openmc_particle_restart
!===============================================================================
! READ_PARTICLE_RESTART reads the particle restart file

70
src/plot.F90
View file

@ -2,8 +2,6 @@ module plot
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use error, only: fatal_error, write_message
use geometry, only: find_cell, check_cell_overlap
@ -32,7 +30,8 @@ contains
! RUN_PLOT controls the logic for making one or many plots
!===============================================================================
subroutine openmc_plot_geometry() bind(C)
function openmc_plot_geometry() result(err) bind(C)
integer(C_INT) :: err
integer :: i ! loop index for plots
@ -52,7 +51,8 @@ contains
end associate
end do
end subroutine openmc_plot_geometry
err = 0
end function openmc_plot_geometry
!===============================================================================
! POSITION_RGB computes the red/green/blue values for a given plot with the
@ -340,23 +340,46 @@ contains
subroutine create_voxel(pl)
type(ObjectPlot), intent(in) :: pl
integer :: x, y, z ! voxel location indices
integer(C_INT) :: x, y, z ! voxel location indices
integer :: rgb(3) ! colors (red, green, blue) from 0-255
integer :: id ! id of cell or material
integer :: hdf5_err
integer, target :: data(pl%pixels(3),pl%pixels(2))
integer(C_INT), target :: data(pl%pixels(3),pl%pixels(2))
integer(HID_T) :: file_id
integer(HID_T) :: dspace
integer(HID_T) :: memspace
integer(HID_T) :: dset
integer(HSIZE_T) :: dims(3)
integer(HSIZE_T) :: dims_slab(3)
integer(HSIZE_T) :: offset(3)
real(8) :: vox(3) ! x, y, and z voxel widths
real(8) :: ll(3) ! lower left starting point for each sweep direction
type(Particle) :: p
type(ProgressBar) :: progress
type(c_ptr) :: f_ptr
interface
subroutine voxel_init(file_id, dims, dspace, dset, memspace) bind(C)
import HID_T, HSIZE_T
integer(HID_T), value :: file_id
integer(HSIZE_T), intent(in) :: dims(*)
integer(HID_T), intent(out) :: dspace
integer(HID_T), intent(out) :: dset
integer(HID_T), intent(out) :: memspace
end subroutine voxel_init
subroutine voxel_write_slice(x, dspace, dset, memspace, buf) bind(C)
import C_INT, HID_T, C_PTR
integer(C_INT), value :: x
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
type(C_PTR), value :: buf
end subroutine voxel_write_slice
subroutine voxel_finalize(dspace, dset, memspace) bind(C)
import HID_T
integer(HID_T), value :: dspace
integer(HID_T), value :: dset
integer(HID_T), value :: memspace
end subroutine voxel_finalize
end interface
! compute voxel widths in each direction
vox = pl % width/dble(pl % pixels)
@ -371,7 +394,7 @@ contains
p % coord(1) % universe = root_universe
! Open binary plot file for writing
file_id = file_create(pl%path_plot)
file_id = file_open(pl%path_plot, 'w')
! write header info
call write_attribute(file_id, "filetype", 'voxel')
@ -390,20 +413,8 @@ contains
! Create dataset for voxel data -- note that the dimensions are reversed
! since we want the order in the file to be z, y, x
dims(:) = [pl%pixels(3), pl%pixels(2), pl%pixels(1)]
call h5screate_simple_f(3, dims, dspace, hdf5_err)
call h5dcreate_f(file_id, "data", H5T_NATIVE_INTEGER, dspace, dset, hdf5_err)
! Create another dataspace for 2D array in memory
dims_slab(1) = pl%pixels(3)
dims_slab(2) = pl%pixels(2)
dims_slab(3) = 1
call h5screate_simple_f(2, dims_slab(1:2), memspace, hdf5_err)
! Initialize offset and get pointer to data
offset(:) = 0
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims_slab, hdf5_err)
f_ptr = c_loc(data)
dims(:) = pl % pixels
call voxel_init(file_id, dims, dspace, dset, memspace)
! move to center of voxels
ll = ll + vox / TWO
@ -433,15 +444,10 @@ contains
p % coord(1) % xyz(3) = ll(3)
! Write to HDF5 dataset
offset(3) = x - 1
call h5soffset_simple_f(dspace, offset, hdf5_err)
call h5dwrite_f(dset, H5T_NATIVE_INTEGER, f_ptr, hdf5_err, &
mem_space_id=memspace, file_space_id=dspace)
call voxel_write_slice(x, dspace, dset, memspace, c_loc(data))
end do
call h5dclose_f(dset, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call voxel_finalize(dspace, dset, memspace)
call file_close(file_id)
end subroutine create_voxel

42
src/plot.cpp Normal file
View file

@ -0,0 +1,42 @@
#include "plot.h"
namespace openmc {
void
voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace, hid_t* dset,
hid_t* memspace)
{
// Create dataspace/dataset for voxel data
*dspace = H5Screate_simple(3, dims, nullptr);
*dset = H5Dcreate(file_id, "data", H5T_NATIVE_INT, *dspace, H5P_DEFAULT,
H5P_DEFAULT, H5P_DEFAULT);
// Create dataspace for a slice of the voxel
hsize_t dims_slice[2] {dims[1], dims[2]};
*memspace = H5Screate_simple(2, dims_slice, nullptr);
// Select hyperslab in dataspace
hsize_t start[3] {0, 0, 0};
hsize_t count[3] {1, dims[1], dims[2]};
H5Sselect_hyperslab(*dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
}
void
voxel_write_slice(int x, hid_t dspace, hid_t dset, hid_t memspace, void* buf)
{
hssize_t offset[3] {x - 1, 0, 0};
H5Soffset_simple(dspace, offset);
H5Dwrite(dset, H5T_NATIVE_INT, memspace, dspace, H5P_DEFAULT, buf);
}
void
voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
{
H5Dclose(dset);
H5Sclose(dspace);
H5Sclose(memspace);
}
} // namespace openmc

15
src/plot.h Normal file
View file

@ -0,0 +1,15 @@
#ifndef PLOT_H
#define PLOT_H
#include "hdf5.h"
namespace openmc {
extern "C" void voxel_init(hid_t file_id, const hsize_t* dims, hid_t* dspace,
hid_t* dset, hid_t* memspace);
extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
hid_t memspace, void* buf);
extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
} // namespace openmc
#endif // PLOT_H

4
src/product_header.F90
View file

@ -1,13 +1,11 @@
module product_header
use hdf5, only: HID_T
use angleenergy_header, only: AngleEnergyContainer
use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, &
EMISSION_TOTAL, NEUTRON, PHOTON
use endf_header, only: Tabulated1D, Function1D, Polynomial
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, close_dataset, read_dataset
open_dataset, close_dataset, read_dataset, HID_T
use random_lcg, only: prn
use secondary_correlated, only: CorrelatedAngleEnergy
use secondary_kalbach, only: KalbachMann

6
src/reaction_header.F90
View file

@ -1,11 +1,7 @@
module reaction_header
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5gget_info_f, h5lget_name_by_idx_f, &
H5_INDEX_NAME_F, H5_ITER_INC_F
use constants, only: MAX_WORD_LEN
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape, get_groups
use hdf5_interface
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with

9
src/sab_header.F90
View file

@ -8,10 +8,7 @@ module sab_header
use dict_header, only: DictIntInt, DictCharInt
use distribution_univariate, only: Tabular
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets, object_exists, &
get_name
use hdf5_interface
use random_lcg, only: prn
use secondary_correlated, only: CorrelatedAngleEnergy
use settings
@ -103,7 +100,6 @@ contains
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer(SIZE_T) :: name_len
integer(HID_T) :: T_group
integer(HID_T) :: elastic_group
integer(HID_T) :: inelastic_group
@ -123,8 +119,7 @@ contains
type(VectorInt) :: temps_to_read
! Get name of table from group
name_len = len(this % name)
this % name = get_name(group_id, name_len)
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))

5
src/secondary_correlated.F90
View file

@ -1,13 +1,10 @@
module secondary_correlated
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, ONE, HALF, TWO, HISTOGRAM, LINEAR_LINEAR
use distribution_univariate, only: DistributionContainer, Tabular
use hdf5_interface, only: get_shape, read_attribute, open_dataset, &
read_dataset, close_dataset
use hdf5_interface
use random_lcg, only: prn
!===============================================================================

5
src/secondary_kalbach.F90
View file

@ -1,12 +1,9 @@
module secondary_kalbach
use hdf5, only: HID_T, HSIZE_T
use algorithm, only: binary_search
use angleenergy_header, only: AngleEnergy
use constants, only: ZERO, HALF, ONE, TWO, HISTOGRAM, LINEAR_LINEAR
use hdf5_interface, only: read_attribute, read_dataset, open_dataset, &
close_dataset, get_shape
use hdf5_interface
use random_lcg, only: prn
!===============================================================================

4
src/secondary_nbody.F90
View file

@ -1,10 +1,8 @@
module secondary_nbody
use hdf5, only: HID_T
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, PI
use hdf5_interface, only: read_attribute
use hdf5_interface, only: read_attribute, HID_T
use math, only: maxwell_spectrum
use random_lcg, only: prn

4
src/secondary_uncorrelated.F90
View file

@ -1,7 +1,5 @@
module secondary_uncorrelated
use hdf5, only: HID_T
use angle_distribution, only: AngleDistribution
use angleenergy_header, only: AngleEnergy
use constants, only: ONE, TWO, MAX_WORD_LEN
@ -9,7 +7,7 @@ module secondary_uncorrelated
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton
use error, only: warning
use hdf5_interface, only: read_attribute, open_group, close_group, &
object_exists
object_exists, HID_T
use random_lcg, only: prn
!===============================================================================

14
src/settings.F90
View file

@ -75,16 +75,16 @@ module settings
real(8) :: weight_survive = ONE
! Mode to run in (fixed source, eigenvalue, plotting, etc)
integer(C_INT), bind(C) :: run_mode = NONE
integer(C_INT), bind(C, name='openmc_run_mode') :: run_mode = NONE
! Restart run
logical :: restart_run = .false.
logical(C_BOOL), bind(C, name='openmc_restart_run') :: restart_run = .false.
! The verbosity controls how much information will be printed to the screen
! and in logs
integer(C_INT), bind(C) :: verbosity = 7
integer(C_INT), bind(C, name='openmc_verbosity') :: verbosity = 7
logical :: check_overlaps = .false.
logical(C_BOOL), bind(C, name='openmc_check_overlaps') :: check_overlaps = .false.
! Trace for single particle
integer :: trace_batch
@ -92,11 +92,13 @@ module settings
integer(8) :: trace_particle
! Particle tracks
logical :: write_all_tracks = .false.
logical(C_BOOL), bind(C, name='openmc_write_all_tracks') :: &
write_all_tracks = .false.
integer, allocatable :: track_identifiers(:,:)
! Particle restart run
logical :: particle_restart_run = .false.
logical(C_BOOL), bind(C, name='openmc_particle_restart_run') :: &
particle_restart_run = .false.
! Write out initial source
logical :: write_initial_source = .false.

64
src/simulation.F90
View file

@ -44,40 +44,32 @@ module simulation
implicit none
private
public :: openmc_next_batch
public :: openmc_run
public :: openmc_simulation_init
public :: openmc_simulation_finalize
integer(C_INT), parameter :: STATUS_EXIT_NORMAL = 0
integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1
integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2
contains
!===============================================================================
! OPENMC_RUN encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a fixed source or k-eigenvalue
! calculation.
!===============================================================================
subroutine openmc_run() bind(C)
call openmc_simulation_init()
do while (openmc_next_batch() == 0)
end do
call openmc_simulation_finalize()
end subroutine openmc_run
!===============================================================================
! OPENMC_NEXT_BATCH
!===============================================================================
function openmc_next_batch() result(retval) bind(C)
integer(C_INT) :: retval
function openmc_next_batch(status) result(err) bind(C)
integer(C_INT), intent(out), optional :: status
integer(C_INT) :: err
type(Particle) :: p
integer(8) :: i_work
err = 0
! Make sure simulation has been initialized
if (.not. simulation_initialized) then
retval = -3
err = E_ALLOCATE
call set_errmsg("Simulation has not been initialized yet.")
return
end if
@ -86,7 +78,7 @@ contains
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
retval = 0
status = STATUS_EXIT_NORMAL
return
end if
@ -124,12 +116,14 @@ contains
call finalize_batch()
! Check simulation ending criteria
if (current_batch == n_max_batches) then
retval = -1
elseif (satisfy_triggers) then
retval = -2
else
retval = 0
if (present(status)) then
if (current_batch == n_max_batches) then
status = STATUS_EXIT_MAX_BATCH
elseif (satisfy_triggers) then
status = STATUS_EXIT_ON_TRIGGER
else
status = STATUS_EXIT_NORMAL
end if
end if
end function openmc_next_batch
@ -318,6 +312,7 @@ contains
subroutine finalize_batch()
integer(C_INT) :: err
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
@ -353,7 +348,7 @@ contains
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call openmc_statepoint_write()
err = openmc_statepoint_write()
end if
! Write out source point if it's been specified for this batch
@ -402,9 +397,13 @@ contains
! INITIALIZE_SIMULATION
!===============================================================================
subroutine openmc_simulation_init() bind(C)
function openmc_simulation_init() result(err) bind(C)
integer(C_INT) :: err
integer :: i
err = 0
! Skip if simulation has already been initialized
if (simulation_initialized) return
@ -462,14 +461,15 @@ contains
! Set flag indicating initialization is done
simulation_initialized = .true.
end subroutine openmc_simulation_init
end function openmc_simulation_init
!===============================================================================
! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays
! execution time and results
!===============================================================================
subroutine openmc_simulation_finalize() bind(C)
function openmc_simulation_finalize() result(err) bind(C)
integer(C_INT) :: err
integer :: i ! loop index
#ifdef OPENMC_MPI
@ -485,6 +485,8 @@ contains
#endif
#endif
err = 0
! Skip if simulation was never run
if (.not. simulation_initialized) return
@ -559,7 +561,7 @@ contains
need_depletion_rx = .false.
simulation_initialized = .false.
end subroutine openmc_simulation_finalize
end function openmc_simulation_finalize
!===============================================================================
! CALCULATE_WORK determines how many particles each processor should simulate

18
src/simulation.cpp Normal file
View file

@ -0,0 +1,18 @@
#include "openmc.h"
// OPENMC_RUN encompasses all the main logic where iterations are performed
// over the batches, generations, and histories in a fixed source or k-eigenvalue
// calculation.
int openmc_run() {
openmc_simulation_init();
int err = 0;
int status = 0;
while (status == 0 && err == 0) {
err = openmc_next_batch(&status);
}
openmc_simulation_finalize();
return err;
}

13
src/simulation_header.F90
View file

@ -23,7 +23,8 @@ module simulation_header
integer :: current_batch ! current batch
integer :: current_gen ! current generation within a batch
integer :: total_gen = 0 ! total number of generations simulated
logical(C_BOOL), bind(C) :: simulation_initialized = .false.
logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: &
simulation_initialized = .false.
logical :: need_depletion_rx ! need to calculate depletion reaction rx?
! ============================================================================
@ -31,8 +32,8 @@ module simulation_header
logical :: satisfy_triggers = .false. ! whether triggers are satisfied
integer(8) :: work ! number of particles per processor
integer(8), allocatable :: work_index(:) ! starting index in source bank for each process
integer(C_INT64_T), bind(C, name='openmc_work') :: work ! number of particles per processor
integer(C_INT64_T), allocatable :: work_index(:) ! starting index in source bank for each process
integer(8) :: current_work ! index in source bank of current history simulated
! ============================================================================
@ -40,8 +41,8 @@ module simulation_header
! Temporary k-effective values
type(VectorReal) :: k_generation ! single-generation estimates of k
real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches
real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k
real(C_DOUBLE), bind(C, name='openmc_keff') :: keff = ONE ! average k over active batches
real(C_DOUBLE), bind(C, name='openmc_keff_std') :: keff_std ! standard deviation of average k
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
@ -57,7 +58,7 @@ module simulation_header
! PARALLEL PROCESSING VARIABLES
#ifdef _OPENMP
integer :: n_threads = NONE ! number of OpenMP threads
integer(C_INT), bind(C, name='openmc_n_threads') :: n_threads = NONE ! number of OpenMP threads
integer :: thread_id ! ID of a given thread
#endif

13
src/source.F90
View file

@ -1,6 +1,5 @@
module source
use hdf5, only: HID_T
#ifdef OPENMC_MPI
use message_passing
#endif
@ -12,7 +11,7 @@ module source
use distribution_multivariate, only: SpatialBox
use error, only: fatal_error
use geometry, only: find_cell
use hdf5_interface, only: file_create, file_open, file_close, read_dataset
use hdf5_interface
use math
use message_passing, only: rank
use mgxs_header, only: rev_energy_bins, num_energy_groups
@ -55,16 +54,16 @@ contains
file_id = file_open(path_source, 'r', parallel=.true.)
! Read the file type
call read_dataset(filetype, file_id, "filetype")
call read_attribute(filetype, file_id, "filetype")
! Check to make sure this is a source file
if (filetype /= 'source') then
if (filetype /= 'source' .and. filetype /= 'statepoint') then
call fatal_error("Specified starting source file not a source file &
&type.")
end if
! Read in the source bank
call read_source_bank(file_id)
call read_source_bank(file_id, work_index, source_bank)
! Close file
call file_close(file_id)
@ -88,8 +87,8 @@ contains
if (write_initial_source) then
call write_message('Writing out initial source...', 5)
filename = trim(path_output) // 'initial_source.h5'
file_id = file_create(filename, parallel=.true.)
call write_source_bank(file_id)
file_id = file_open(filename, 'w', parallel=.true.)
call write_source_bank(file_id, work_index, source_bank)
call file_close(file_id)
end if

266
src/state_point.F90
View file

@ -13,8 +13,7 @@ module state_point
use, intrinsic :: ISO_C_BINDING
use hdf5
use bank_header, only: Bank
use cmfd_header
use constants
use eigenvalue, only: openmc_get_keff
@ -37,30 +36,45 @@ module state_point
implicit none
interface
subroutine write_source_bank(group_id, work_index, bank_) bind(C)
import HID_T, C_INT64_T, Bank
integer(HID_T), value :: group_id
integer(C_INT64_T), intent(in) :: work_index(*)
type(Bank), intent(in) :: bank_(*)
end subroutine write_source_bank
subroutine read_source_bank(group_id, work_index, bank_) bind(C)
import HID_T, C_INT64_T, Bank
integer(HID_T), value :: group_id
integer(C_INT64_T), intent(in) :: work_index(*)
type(Bank), intent(out) :: bank_(*)
end subroutine read_source_bank
end interface
contains
!===============================================================================
! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk
!===============================================================================
subroutine openmc_statepoint_write(filename) bind(C)
function openmc_statepoint_write(filename) result(err) bind(C)
type(C_PTR), intent(in), optional :: filename
integer(C_INT) :: err
integer :: i, j, k
integer :: i_xs
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer, allocatable :: id_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
integer(C_INT) :: err
real(C_DOUBLE) :: k_combined(2)
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(C_CHAR), pointer :: string(:)
character(len=:, kind=C_CHAR), allocatable :: filename_
err = 0
if (present(filename)) then
call c_f_pointer(filename, string, [MAX_FILE_LEN])
filename_ = to_f_string(string)
@ -76,7 +90,7 @@ contains
if (master) then
! Create statepoint file
file_id = file_create(filename_)
file_id = file_open(filename_, 'w')
! Write file type
call write_attribute(file_id, "filetype", "statepoint")
@ -320,42 +334,9 @@ contains
str_array(j) = reaction_name(tally % score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
call write_dataset(tally_group, "n_user_score_bins", &
tally % n_user_score_bins)
deallocate(str_array)
! Write explicit moment order strings for each score bin
k = 1
allocate(str_array(tally % n_score_bins))
MOMENT_LOOP: do j = 1, tally % n_user_score_bins
select case(tally % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
str_array(k) = trim(to_str(tally % moment_order(k)))
k = k + 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
do n_order = 0, tally % moment_order(k)
str_array(k) = trim(to_str(n_order))
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
str_array(k) = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
k = k + 1
end do
end do
case default
str_array(k) = ''
k = k + 1
end select
end do MOMENT_LOOP
call write_dataset(tally_group, "moment_orders", str_array)
deallocate(str_array)
call close_group(tally_group)
end associate
end do TALLY_METADATA
@ -443,7 +424,7 @@ contains
call file_close(file_id)
end if
end subroutine openmc_statepoint_write
end function openmc_statepoint_write
!===============================================================================
! WRITE_SOURCE_POINT
@ -476,8 +457,8 @@ contains
! Create separate source file
if (master .or. parallel) then
file_id = file_create(filename, parallel=.true.)
call write_dataset(file_id, "filetype", 'source')
file_id = file_open(filename, 'w', parallel=.true.)
call write_attribute(file_id, "filetype", 'source')
end if
else
filename = trim(path_output) // 'statepoint.' // &
@ -485,11 +466,11 @@ contains
filename = trim(filename) // '.h5'
if (master .or. parallel) then
file_id = file_open(filename, 'w', parallel=.true.)
file_id = file_open(filename, 'a', parallel=.true.)
end if
end if
call write_source_bank(file_id)
call write_source_bank(file_id, work_index, source_bank)
if (master .or. parallel) call file_close(file_id)
end if
@ -498,11 +479,11 @@ contains
filename = trim(path_output) // 'source' // '.h5'
call write_message("Creating source file " // trim(filename) // "...", 5)
if (master .or. parallel) then
file_id = file_create(filename, parallel=.true.)
call write_dataset(file_id, "filetype", 'source')
file_id = file_open(filename, 'w', parallel=.true.)
call write_attribute(file_id, "filetype", 'source')
end if
call write_source_bank(file_id)
call write_source_bank(file_id, work_index, source_bank)
if (master .or. parallel) call file_close(file_id)
end if
@ -815,14 +796,10 @@ contains
! Open source file
file_id = file_open(path_source_point, 'r', parallel=.true.)
! Read file type
call read_dataset(int_array(1), file_id, "filetype")
end if
! Write out source
call read_source_bank(file_id)
call read_source_bank(file_id, work_index, source_bank)
end if
@ -831,187 +808,4 @@ contains
end subroutine load_state_point
!===============================================================================
! WRITE_SOURCE_BANK writes OpenMC source_bank data
!===============================================================================
subroutine write_source_bank(group_id)
use bank_header, only: Bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data or file space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: offset(1) ! source data offset
integer(HSIZE_T) :: dims(1)
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#else
integer :: i
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
type(Bank), allocatable, target :: temp_source(:)
#endif
#endif
#ifdef PHDF5
! Set size of total dataspace for all procs and rank
dims(1) = n_particles
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, dspace, dset, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(1, dims, memspace, hdf5_err)
! Select hyperslab for this dataspace
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up the property list for parallel writing
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
! Write data to file in parallel
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
call h5pclose_f(plist, hdf5_err)
#else
if (master) then
! Create dataset big enough to hold all source sites
dims(1) = n_particles
call h5screate_simple_f(1, dims, dspace, hdf5_err)
call h5dcreate_f(group_id, "source_bank", hdf5_bank_t, &
dspace, dset, hdf5_err)
! Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
allocate(temp_source(work))
temp_source(:) = source_bank(:)
#endif
do i = 0, n_procs - 1
! Create memory space
dims(1) = work_index(i+1) - work_index(i)
call h5screate_simple_f(1, dims, memspace, hdf5_err)
#ifdef OPENMC_MPI
! Receive source sites from other processes
if (i > 0) then
call MPI_RECV(source_bank, int(dims(1)), MPI_BANK, i, i, &
mpi_intracomm, MPI_STATUS_IGNORE, mpi_err)
end if
#endif
! Select hyperslab for this dataspace
call h5dget_space_f(dset, dspace, hdf5_err)
offset(1) = work_index(i)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up pointer to data and write data to hyperslab
f_ptr = c_loc(source_bank)
call h5dwrite_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace)
call h5sclose_f(memspace, hdf5_err)
call h5sclose_f(dspace, hdf5_err)
end do
! Close all ids
call h5dclose_f(dset, hdf5_err)
! Restore state of source bank
#ifdef OPENMC_MPI
source_bank(:) = temp_source(:)
deallocate(temp_source)
#endif
else
#ifdef OPENMC_MPI
call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, &
mpi_intracomm, mpi_err)
#endif
end if
#endif
end subroutine write_source_bank
!===============================================================================
! READ_SOURCE_BANK reads OpenMC source_bank data
!===============================================================================
subroutine read_source_bank(group_id)
use bank_header, only: Bank
integer(HID_T), intent(in) :: group_id
integer :: hdf5_err
integer(HID_T) :: dset ! data set handle
integer(HID_T) :: dspace ! data space handle
integer(HID_T) :: memspace ! memory space handle
integer(HSIZE_T) :: dims(1) ! dimensions on one processor
integer(HSIZE_T) :: dims_all(1) ! overall dimensions
integer(HSIZE_T) :: maxdims(1) ! maximum dimensions
integer(HSIZE_T) :: offset(1) ! offset of data
type(c_ptr) :: f_ptr
#ifdef PHDF5
integer(HID_T) :: plist ! property list
#endif
! Open the dataset
call h5dopen_f(group_id, "source_bank", dset, hdf5_err)
! Create another data space but for each proc individually
dims(1) = work
call h5screate_simple_f(1, dims, memspace, hdf5_err)
! Make sure source bank is big enough
call h5dget_space_f(dset, dspace, hdf5_err)
call h5sget_simple_extent_dims_f(dspace, dims_all, maxdims, hdf5_err)
if (size(source_bank, KIND=HSIZE_T) > dims_all(1)) then
call fatal_error("Number of source sites in source file is less than &
&number of source particles per generation.")
end if
! Select hyperslab for each process
offset(1) = work_index(rank)
call h5sselect_hyperslab_f(dspace, H5S_SELECT_SET_F, offset, dims, hdf5_err)
! Set up pointer to data
f_ptr = c_loc(source_bank)
#ifdef PHDF5
! Read data in parallel
call h5pcreate_f(H5P_DATASET_XFER_F, plist, hdf5_err)
call h5pset_dxpl_mpio_f(plist, H5FD_MPIO_COLLECTIVE_F, hdf5_err)
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace, &
xfer_prp=plist)
call h5pclose_f(plist, hdf5_err)
#else
call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err, &
file_space_id=dspace, mem_space_id=memspace)
#endif
! Close all ids
call h5sclose_f(dspace, hdf5_err)
call h5sclose_f(memspace, hdf5_err)
call h5dclose_f(dset, hdf5_err)
end subroutine read_source_bank
end module state_point

165
src/state_point.cpp Normal file
View file

@ -0,0 +1,165 @@
#include "state_point.h"
#include <algorithm>
#include <vector>
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "error.h"
#include "message_passing.h"
#include "openmc.h"
namespace openmc {
hid_t h5banktype() {
// Create type for array of 3 reals
hsize_t dims[] {3};
hid_t triplet = H5Tarray_create(H5T_NATIVE_DOUBLE, 1, dims);
// Create bank datatype
hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Bank));
H5Tinsert(banktype, "wgt", HOFFSET(Bank, wgt), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "xyz", HOFFSET(Bank, xyz), triplet);
H5Tinsert(banktype, "uvw", HOFFSET(Bank, uvw), triplet);
H5Tinsert(banktype, "E", HOFFSET(Bank, E), H5T_NATIVE_DOUBLE);
H5Tinsert(banktype, "delayed_group", HOFFSET(Bank, delayed_group), H5T_NATIVE_INT);
H5Tclose(triplet);
return banktype;
}
void
write_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
{
hid_t banktype = h5banktype();
#ifdef PHDF5
// Set size of total dataspace for all procs and rank
hsize_t dims[] {static_cast<hsize_t>(n_particles)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t count[] {static_cast<hsize_t>(openmc_work)};
hid_t memspace = H5Screate_simple(1, count, nullptr);
// Select hyperslab for this dataspace
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Set up the property list for parallel writing
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
// Write data to file in parallel
H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank);
// Free resources
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Pclose(plist);
#else
if (openmc_master) {
// Create dataset big enough to hold all source sites
hsize_t dims[] {static_cast<hsize_t>(n_particles)};
hid_t dspace = H5Screate_simple(1, dims, nullptr);
hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
// Save source bank sites since the souce_bank array is overwritten below
#ifdef OPENMC_MPI
std::vector<Bank> temp_source {source_bank, source_bank + openmc_work};
#endif
for (int i = 0; i < openmc::mpi::n_procs; ++i) {
// Create memory space
hsize_t count[] {static_cast<hsize_t>(work_index[i+1] - work_index[i])};
hid_t memspace = H5Screate_simple(1, count, nullptr);
#ifdef OPENMC_MPI
// Receive source sites from other processes
if (i > 0)
MPI_Recv(source_bank, count[0], openmc::mpi::bank, i, i,
openmc::mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// Select hyperslab for this dataspace
dspace = H5Dget_space(dset);
hsize_t start[] {static_cast<hsize_t>(work_index[i])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Write data to hyperslab
H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
H5Sclose(memspace);
H5Sclose(dspace);
}
// Close all ids
H5Dclose(dset);
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank);
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank, openmc_work, openmc::mpi::bank, 0, openmc::mpi::rank,
openmc::mpi::intracomm);
#endif
}
#endif
H5Tclose(banktype);
}
void read_source_bank(hid_t group_id, int64_t* work_index, Bank* source_bank)
{
hid_t banktype = h5banktype();
// Open the dataset
hid_t dset = H5Dopen(group_id, "source_bank", H5P_DEFAULT);
// Create another data space but for each proc individually
hsize_t dims[] {static_cast<hsize_t>(openmc_work)};
hid_t memspace = H5Screate_simple(1, dims, nullptr);
// Make sure source bank is big enough
hid_t dspace = H5Dget_space(dset);
hsize_t dims_all[1];
H5Sget_simple_extent_dims(dspace, dims_all, nullptr);
if (work_index[openmc::mpi::n_procs] > dims_all[0]) {
fatal_error("Number of source sites in source file is less "
"than number of source particles per generation.");
}
// Select hyperslab for each process
hsize_t start[] {static_cast<hsize_t>(work_index[openmc::mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, dims, nullptr);
#ifdef PHDF5
// Read data in parallel
hid_t plist = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
H5Pclose(plist);
#else
H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
#endif
// Close all ids
H5Sclose(dspace);
H5Sclose(memspace);
H5Dclose(dset);
H5Tclose(banktype);
}
} // namespace openmc

17
src/state_point.h Normal file
View file

@ -0,0 +1,17 @@
#ifndef STATE_POINT_H
#define STATE_POINT_H
#include <cstdint>
#include "hdf5.h"
#include "openmc.h"
namespace openmc {
extern "C" void write_source_bank(hid_t group_id, int64_t* work_index,
Bank* source_bank);
extern "C" void read_source_bank(hid_t group_id, int64_t* work_index,
Bank* source_bank);
} // namespace openmc
#endif // STATE_POINT_H

158
src/summary.F90
View file

@ -1,7 +1,5 @@
module summary
use hdf5
use constants
use endf, only: reaction_name
use error, only: write_message
@ -38,7 +36,7 @@ contains
call write_message("Writing summary.h5 file...", 5)
! Create a new file using default properties.
file_id = file_create("summary.h5")
file_id = file_open("summary.h5", 'w')
call write_header(file_id)
call write_nuclides(file_id)
@ -77,34 +75,81 @@ contains
subroutine write_nuclides(file_id)
integer(HID_T), intent(in) :: file_id
integer(HID_T) :: nuclide_group
integer(HID_T) :: macro_group
integer :: i
character(12), allocatable :: nucnames(:)
character(12), allocatable :: nuc_names(:)
character(12), allocatable :: macro_names(:)
real(8), allocatable :: awrs(:)
integer :: num_nuclides
integer :: num_macros
integer :: j
integer :: k
! Write useful data from nuclide objects
nuclide_group = create_group(file_id, "nuclides")
call write_attribute(nuclide_group, "n_nuclides", n_nuclides)
! Find how many of these nuclides are macroscopic objects
if (run_CE) then
! Then none are macroscopic
num_nuclides = n_nuclides
num_macros = 0
else
num_nuclides = 0
num_macros = 0
do i = 1, n_nuclides
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
end if
end do
end if
! Build array of nuclide names and awrs
allocate(nucnames(n_nuclides))
allocate(awrs(n_nuclides))
! Build array of nuclide names and awrs while only sorting nuclides from
! macroscopics
if (num_nuclides > 0) then
allocate(nuc_names(num_nuclides))
allocate(awrs(num_nuclides))
end if
if (num_macros > 0) then
allocate(macro_names(num_macros))
end if
j = 1
k = 1
do i = 1, n_nuclides
if (run_CE) then
nucnames(i) = nuclides(i) % name
nuc_names(i) = nuclides(i) % name
awrs(i) = nuclides(i) % awr
else
nucnames(i) = nuclides_MG(i) % obj % name
awrs(i) = nuclides_MG(i) % obj % awr
if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(j) = nuclides_MG(i) % obj % name
awrs(j) = nuclides_MG(i) % obj % awr
j = j + 1
else
macro_names(k) = nuclides_MG(i) % obj % name
k = k + 1
end if
end if
end do
nuclide_group = create_group(file_id, "nuclides")
call write_attribute(nuclide_group, "n_nuclides", num_nuclides)
macro_group = create_group(file_id, "macroscopics")
call write_attribute(macro_group, "n_macroscopics", num_macros)
! Write nuclide names and awrs
call write_dataset(nuclide_group, "names", nucnames)
call write_dataset(nuclide_group, "awrs", awrs)
if (num_nuclides > 0) then
! Write useful data from nuclide objects
call write_dataset(nuclide_group, "names", nuc_names)
call write_dataset(nuclide_group, "awrs", awrs)
end if
if (num_macros > 0) then
! Write useful data from macroscopic objects
call write_dataset(macro_group, "names", macro_names)
end if
call close_group(nuclide_group)
call close_group(macro_group)
deallocate(nucnames, awrs)
if (allocated(nuc_names)) deallocate(nuc_names, awrs)
if (allocated(macro_names)) deallocate(macro_names)
end subroutine write_nuclides
@ -371,7 +416,13 @@ contains
integer :: i
integer :: j
character(20), allocatable :: nucnames(:)
integer :: k
integer :: n
character(20), allocatable :: nuc_names(:)
character(20), allocatable :: macro_names(:)
real(8), allocatable :: nuc_densities(:)
integer :: num_nuclides
integer :: num_macros
integer(HID_T) :: materials_group
integer(HID_T) :: material_group
type(Material), pointer :: m
@ -399,24 +450,67 @@ contains
! Write atom density with units
call write_dataset(material_group, "atom_density", m % density)
! Copy ZAID for each nuclide to temporary array
allocate(nucnames(m%n_nuclides))
do j = 1, m%n_nuclides
if (run_CE) then
nucnames(j) = nuclides(m%nuclide(j))%name
else
nucnames(j) = nuclides_MG(m%nuclide(j))%obj%name
if (run_CE) then
num_nuclides = m % n_nuclides
num_macros = 0
else
! Find the number of macroscopic and nuclide data in this material
num_nuclides = 0
num_macros = 0
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
num_nuclides = num_nuclides + 1
else
num_macros = num_macros + 1
end if
end do
end if
! Copy ZAID or macro name for each nuclide to temporary array
if (num_nuclides > 0) then
allocate(nuc_names(num_nuclides))
allocate(nuc_densities(num_nuclides))
end if
if (run_CE) then
do j = 1, m % n_nuclides
nuc_names(j) = nuclides(m%nuclide(j))%name
nuc_densities(j) = m % atom_density(j)
end do
else
if (num_macros > 0) then
allocate(macro_names(num_macros))
end if
end do
k = 1
n = 1
do j = 1, m % n_nuclides
if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then
nuc_names(k) = nuclides_MG(m % nuclide(j)) % obj % name
nuc_densities(k) = m % atom_density(j)
k = k + 1
else
macro_names(n) = nuclides_MG(m % nuclide(j)) % obj % name
n = n + 1
end if
end do
end if
! Write temporary array to 'nuclides'
call write_dataset(material_group, "nuclides", nucnames)
if (num_nuclides > 0) then
call write_dataset(material_group, "nuclides", nuc_names)
! Deallocate temporary array
deallocate(nuc_names)
! Write atom densities
call write_dataset(material_group, "nuclide_densities", nuc_densities)
deallocate(nuc_densities)
end if
! Deallocate temporary array
deallocate(nucnames)
! Write atom densities
call write_dataset(material_group, "nuclide_densities", m%atom_density)
! Write temporary array to 'macroscopics'
if (num_macros > 0) then
call write_dataset(material_group, "macroscopics", macro_names)
! Deallocate temporary array
deallocate(macro_names)
end if
if (m%n_sab > 0) then
call write_dataset(material_group, "sab_names", m%sab_names)

83
src/surface.cpp
View file

@ -11,6 +11,12 @@
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
int32_t n_surfaces;
//==============================================================================
// Helper functions for reading the "coeffs" node of an XML surface element
//==============================================================================
@ -202,21 +208,21 @@ Surface::to_hdf5(hid_t group_id) const
switch(bc) {
case BC_TRANSMIT :
write_string(surf_group, "boundary_type", "transmission");
write_string(surf_group, "boundary_type", "transmission", false);
break;
case BC_VACUUM :
write_string(surf_group, "boundary_type", "vacuum");
write_string(surf_group, "boundary_type", "vacuum", false);
break;
case BC_REFLECT :
write_string(surf_group, "boundary_type", "reflective");
write_string(surf_group, "boundary_type", "reflective", false);
break;
case BC_PERIODIC :
write_string(surf_group, "boundary_type", "periodic");
write_string(surf_group, "boundary_type", "periodic", false);
break;
}
if (!name.empty()) {
write_string(surf_group, "name", name);
write_string(surf_group, "name", name, false);
}
to_hdf5_inner(surf_group);
@ -297,7 +303,7 @@ inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-plane");
write_string(group_id, "type", "x-plane", false);
std::array<double, 1> coeffs {{x0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -362,7 +368,7 @@ inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-plane");
write_string(group_id, "type", "y-plane", false);
std::array<double, 1> coeffs {{y0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -428,7 +434,7 @@ inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const
void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-plane");
write_string(group_id, "type", "z-plane", false);
std::array<double, 1> coeffs {{z0}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -489,7 +495,7 @@ SurfacePlane::normal(const double xyz[3], double uvw[3]) const
void SurfacePlane::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "plane");
write_string(group_id, "type", "plane", false);
std::array<double, 4> coeffs {{A, B, C, D}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -621,7 +627,7 @@ inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cylinder");
write_string(group_id, "type", "x-cylinder", false);
std::array<double, 3> coeffs {{y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -655,7 +661,7 @@ inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cylinder");
write_string(group_id, "type", "y-cylinder", false);
std::array<double, 3> coeffs {{x0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -689,7 +695,7 @@ inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const
void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cylinder");
write_string(group_id, "type", "z-cylinder", false);
std::array<double, 3> coeffs {{x0, y0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -760,7 +766,7 @@ inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const
void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "sphere");
write_string(group_id, "type", "sphere", false);
std::array<double, 4> coeffs {{x0, y0, z0, r}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -877,7 +883,7 @@ inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "x-cone");
write_string(group_id, "type", "x-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -911,7 +917,7 @@ inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "y-cone");
write_string(group_id, "type", "y-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -945,7 +951,7 @@ inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const
void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "z-cone");
write_string(group_id, "type", "z-cone", false);
std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -1039,7 +1045,7 @@ SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const
void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
{
write_string(group_id, "type", "quadric");
write_string(group_id, "type", "quadric", false);
std::array<double, 10> coeffs {{A, B, C, D, E, F, G, H, J, K}};
write_double_1D(group_id, "coefficients", coeffs);
}
@ -1232,4 +1238,47 @@ read_surfaces(pugi::xml_node *node)
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
extern "C" int surface_id(Surface *surf) {return surf->id;}
extern "C" int surface_bc(Surface *surf) {return surf->bc;}
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
{return surf->sense(xyz, uvw);}
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
extern "C" double
surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
{return surf->distance(xyz, uvw, coincident);}
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
{surf->to_hdf5(group);}
extern "C" int surface_i_periodic(PeriodicSurface *surf)
{return surf->i_periodic;}
extern "C" bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
extern "C" void free_memory_surfaces_c()
{
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = nullptr;
n_surfaces = 0;
surface_dict.clear();
}
} // namespace openmc

56
src/surface.h
View file

@ -25,15 +25,14 @@ extern "C" const int BC_PERIODIC {3};
//==============================================================================
extern "C" double FP_COINCIDENT;
constexpr double INFTY{std::numeric_limits<double>::max()};
constexpr double INFTY {std::numeric_limits<double>::max()};
constexpr int C_NONE {-1};
//==============================================================================
// Global variables
//==============================================================================
// Braces force n_surfaces to be defined here, not just declared.
extern "C" {int32_t n_surfaces {0};}
extern "C" int32_t n_surfaces;
class Surface;
Surface **surfaces_c;
@ -384,44 +383,19 @@ public:
// Fortran compatibility functions
//==============================================================================
extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
extern "C" int surface_id(Surface *surf) {return surf->id;}
extern "C" int surface_bc(Surface *surf) {return surf->bc;}
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
{return surf->sense(xyz, uvw);}
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
extern "C" double
surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
{return surf->distance(xyz, uvw, coincident);}
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
{surf->to_hdf5(group);}
extern "C" int surface_i_periodic(PeriodicSurface *surf)
{return surf->i_periodic;}
extern "C" bool
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
extern "C" void free_memory_surfaces_c()
{
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = nullptr;
n_surfaces = 0;
surface_dict.clear();
}
extern "C" Surface* surface_pointer(int surf_ind);
extern "C" int surface_id(Surface *surf);
extern "C" int surface_bc(Surface *surf);
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
extern "C" double surface_distance(Surface *surf, double xyz[3], double uvw[3],
bool coincident);
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
extern "C" void surface_to_hdf5(Surface *surf, hid_t group);
extern "C" int surface_i_periodic(PeriodicSurface *surf);
extern "C" bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
double xyz[3], double uvw[3]);
extern "C" void free_memory_surfaces_c();
} // namespace openmc
#endif // SURFACE_H

5
src/surface_header.F90
View file

@ -1,9 +1,9 @@
module surface_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use dict_header, only: DictIntInt
use hdf5_interface
implicit none
@ -71,8 +71,7 @@ module surface_header
subroutine surface_to_hdf5_c(surf_ptr, group) &
bind(C, name='surface_to_hdf5')
use ISO_C_BINDING
use hdf5
import C_PTR, HID_T
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(HID_T), intent(in), value :: group

243
src/tallies/tally.F90
View file

@ -84,7 +84,6 @@ contains
integer :: i ! loop index for scoring bins
integer :: l ! loop index for nuclides in material
integer :: m ! loop index for reactions
integer :: q ! loop index for scoring bins
integer :: i_temp ! temperature index
integer :: i_nuc ! index in nuclides array (from material)
integer :: i_energy ! index in nuclide energy grid
@ -104,9 +103,7 @@ contains
! Pre-collision energy of particle
E = p % last_E
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
i = i + 1
SCORE_LOOP: do i = 1, t % n_score_bins
! determine what type of score bin
score_bin = t % score_bins(i)
@ -120,7 +117,7 @@ contains
select case(score_bin)
case (SCORE_FLUX, SCORE_FLUX_YN)
case (SCORE_FLUX)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to a flux bin. We actually use a collision
! estimator in place of an analog one since there is no way to count
@ -140,7 +137,7 @@ contains
end if
case (SCORE_TOTAL, SCORE_TOTAL_YN)
case (SCORE_TOTAL)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events will score to the total reaction rate. We can just
! use the weight of the particle entering the collision as the
@ -187,7 +184,7 @@ contains
end if
case (SCORE_SCATTER, SCORE_SCATTER_N)
case (SCORE_SCATTER)
if (t % estimator == ESTIMATOR_ANALOG) then
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -197,7 +194,6 @@ contains
score = p % last_wgt * flux
else
! Note SCORE_SCATTER_N not available for tracklength/collision.
if (i_nuclide > 0) then
score = (micro_xs(i_nuclide) % total &
- micro_xs(i_nuclide) % absorption) * atom_density * flux
@ -207,33 +203,7 @@ contains
end if
case (SCORE_SCATTER_PN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + t % moment_order(i)
cycle SCORE_LOOP
end if
! Since only scattering events make it here, again we can use
! the weight entering the collision as the estimator for the
! reaction rate
score = p % last_wgt * flux
case (SCORE_SCATTER_YN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + (t % moment_order(i) + 1)**2 - 1
cycle SCORE_LOOP
end if
! Since only scattering events make it here, again we can use
! the weight entering the collision as the estimator for the
! reaction rate
score = p % last_wgt * flux
case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N)
case (SCORE_NU_SCATTER)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -256,58 +226,6 @@ contains
end if
case (SCORE_NU_SCATTER_PN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + t % moment_order(i)
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p % event_nuclide) % reaction_index(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(E)
end associate
end if
case (SCORE_NU_SCATTER_YN)
! Only analog estimators are available.
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
i = i + (t % moment_order(i) + 1)**2 - 1
cycle SCORE_LOOP
end if
! For scattering production, we need to use the pre-collision
! weight times the yield as the estimate for the number of
! neutrons exiting a reaction with neutrons in the exit channel
if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. &
(p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then
! Don't waste time on very common reactions we know have
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p % event_nuclide) % reaction_index(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(E)
end associate
end if
case (SCORE_ABSORPTION)
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
@ -1281,8 +1199,9 @@ contains
!#########################################################################
! Expand score if necessary and add to tally results.
call expand_and_score(p, t, score_index, filter_index, score_bin, &
score, i)
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
end do SCORE_LOOP
end subroutine score_general_ce
@ -1359,7 +1278,7 @@ contains
end if
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
SCORE_LOOP: do q = 1, t % n_score_bins
i = i + 1
! determine what type of score bin
@ -1374,7 +1293,7 @@ contains
select case(score_bin)
case (SCORE_FLUX, SCORE_FLUX_YN)
case (SCORE_FLUX)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events score to a flux bin. We actually use a collision
! estimator in place of an analog one since there is no way to count
@ -1395,7 +1314,7 @@ contains
end if
case (SCORE_TOTAL, SCORE_TOTAL_YN)
case (SCORE_TOTAL)
if (t % estimator == ESTIMATOR_ANALOG) then
! All events will score to the total reaction rate. We can just
! use the weight of the particle entering the collision as the
@ -1456,15 +1375,10 @@ contains
end if
case (SCORE_SCATTER, SCORE_SCATTER_N, SCORE_SCATTER_PN, SCORE_SCATTER_YN)
case (SCORE_SCATTER)
if (t % estimator == ESTIMATOR_ANALOG) then
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
if (score_bin == SCORE_SCATTER_PN) then
i = i + t % moment_order(i)
else if (score_bin == SCORE_SCATTER_YN) then
i = i + (t % moment_order(i) + 1)**2 - 1
end if
cycle SCORE_LOOP
end if
@ -1485,7 +1399,6 @@ contains
end if
else
! Note SCORE_SCATTER_*N not available for tracklength/collision.
if (i_nuclide > 0) then
score = atom_density * flux * &
nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw)
@ -1498,16 +1411,10 @@ contains
end if
case (SCORE_NU_SCATTER, SCORE_NU_SCATTER_N, SCORE_NU_SCATTER_PN, &
SCORE_NU_SCATTER_YN)
case (SCORE_NU_SCATTER)
if (t % estimator == ESTIMATOR_ANALOG) then
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) then
if (score_bin == SCORE_NU_SCATTER_PN) then
i = i + t % moment_order(i)
else if (score_bin == SCORE_NU_SCATTER_YN) then
i = i + (t % moment_order(i) + 1)**2 - 1
end if
cycle SCORE_LOOP
end if
@ -1528,7 +1435,6 @@ contains
end if
else
! Note SCORE_NU_SCATTER_*N not available for tracklength/collision.
if (i_nuclide > 0) then
score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * &
atom_density * flux
@ -2077,124 +1983,15 @@ contains
!#########################################################################
! Expand score if necessary and add to tally results.
call expand_and_score(p, t, score_index, filter_index, score_bin, &
score, i)
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
end do SCORE_LOOP
nullify(matxs, nucxs)
end subroutine score_general_mg
!===============================================================================
! EXPAND_AND_SCORE takes a previously determined score value and adjusts it
! if necessary (for functional expansion weighting), and then adds the resultant
! value to the tally results array.
!===============================================================================
subroutine expand_and_score(p, t, score_index, filter_index, score_bin, &
score, i)
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(inout) :: score_index
integer, intent(in) :: filter_index ! for % results
integer, intent(in) :: score_bin ! score of concern
real(8), intent(inout) :: score ! data to score
integer, intent(inout) :: i ! Working index
integer :: num_nm ! Number of N,M orders in harmonic
integer :: n ! Moment loop index
real(8) :: uvw(3)
select case(score_bin)
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
! Find the scattering order for a singly requested moment, and
! store its moment contribution.
if (t % moment_order(i) == 1) then
score = score * p % mu ! avoid function call overhead
else
score = score * calc_pn(t % moment_order(i), p % mu)
endif
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
case(SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN)
score_index = score_index - 1
num_nm = 1
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
! multiply score by the angular flux moments and store
!$omp critical (score_general_scatt_yn)
t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, &
filter_index) = t % results(RESULT_VALUE, &
score_index: score_index + num_nm - 1, filter_index) &
+ score * calc_pn(n, p % mu) * calc_rn(n, p % last_uvw)
!$omp end critical (score_general_scatt_yn)
end do
i = i + (t % moment_order(i) + 1)**2 - 1
case(SCORE_FLUX_YN, SCORE_TOTAL_YN)
score_index = score_index - 1
num_nm = 1
if (t % estimator == ESTIMATOR_ANALOG .or. &
t % estimator == ESTIMATOR_COLLISION) then
uvw = p % last_uvw
else if (t % estimator == ESTIMATOR_TRACKLENGTH) then
uvw = p % coord(1) % uvw
end if
! Find the order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + num_nm
! Update number of total n,m bins for this n (m = [-n: n])
num_nm = 2 * n + 1
! multiply score by the angular flux moments and store
!$omp critical (score_general_flux_tot_yn)
t % results(RESULT_VALUE, score_index: score_index + num_nm - 1, &
filter_index) = t % results(RESULT_VALUE, &
score_index: score_index + num_nm - 1, filter_index) &
+ score * calc_rn(n, uvw)
!$omp end critical (score_general_flux_tot_yn)
end do
i = i + (t % moment_order(i) + 1)**2 - 1
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
! Find the scattering order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % moment_order(i)
! determine scoring bin index
score_index = score_index + 1
! get the score and tally it
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) &
+ score * calc_pn(n, p % mu)
end do
i = i + t % moment_order(i)
case default
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
end select
end subroutine expand_and_score
!===============================================================================
! SCORE_ALL_NUCLIDES tallies individual nuclide reaction rates specifically when
! the user requests <nuclides>all</nuclides>.
@ -3134,7 +2931,7 @@ contains
! Currently only one score type
k = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
SCORE_LOOP: do q = 1, t % n_score_bins
k = k + 1
! determine what type of score bin
@ -3144,8 +2941,10 @@ contains
score_index = q
! Expand score if necessary and add to tally results.
call expand_and_score(p, t, score_index, filter_index, score_bin, &
score, k)
!$omp atomic
t % results(RESULT_VALUE, score_index, filter_index) = &
t % results(RESULT_VALUE, score_index, filter_index) + score
end do SCORE_LOOP
! ======================================================================

2
src/tallies/tally_filter.F90
View file

@ -2,8 +2,6 @@ module tally_filter
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use error
use string, only: to_f_string
use tally_filter_header

2
src/tallies/tally_filter_azimuthal.F90
View file

@ -2,8 +2,6 @@ module tally_filter_azimuthal
use, intrinsic :: ISO_C_BINDING
use hdf5
use algorithm, only: binary_search
use constants
use error, only: fatal_error

2
src/tallies/tally_filter_cell.F90
View file

@ -2,8 +2,6 @@ module tally_filter_cell
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error

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