Merge remote-tracking branch 'upstream/develop' into cpp_geometry

This commit is contained in:
Sterling Harper 2018-08-24 20:45:52 -04:00
commit 01b216b63a
145 changed files with 3018 additions and 2325 deletions

View file

@ -17,6 +17,7 @@ cache:
directories:
- $HOME/nndc_hdf5
- $HOME/endf-b-vii.1
- $HOME/WMP_Library
env:
global:
- FC=gfortran
@ -25,14 +26,17 @@ env:
- OMP_NUM_THREADS=2
- OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- OPENMC_ENDF_DATA=$HOME/endf-b-vii.1
- OPENMC_MULTIPOLE_LIBRARY=$HOME/multipole_lib
- OPENMC_MULTIPOLE_LIBRARY=$HOME/WMP_Library
- PATH=$PATH:$HOME/NJOY2016/build
- DISPLAY=:99.0
- COVERALLS_PARALLEL=true
matrix:
- OMP=n MPI=n PHDF5=n
- OMP=y MPI=n PHDF5=n
- OMP=n MPI=y PHDF5=n
- OMP=n MPI=y PHDF5=y
notifications:
webhooks: https://coveralls.io/webhook?repo_token=$COVERALLS_REPO_TOKEN
install:
- ./tools/ci/travis-install.sh
before_script:

View file

@ -288,7 +288,6 @@ set_target_properties(faddeeva PROPERTIES
add_library(libopenmc SHARED
src/algorithm.F90
src/angleenergy_header.F90
src/bank_header.F90
src/api.F90
src/cmfd_data.F90
@ -338,7 +337,6 @@ add_library(libopenmc SHARED
src/reaction_header.F90
src/relaxng
src/sab_header.F90
src/secondary_correlated.F90
src/set_header.F90
src/settings.F90
src/simulation_header.F90
@ -422,11 +420,11 @@ add_library(libopenmc SHARED
src/string_functions.cpp
src/summary.cpp
src/surface.cpp
src/thermal.cpp
src/xml_interface.cpp
src/xsdata.cpp)
set_target_properties(libopenmc PROPERTIES
OUTPUT_NAME openmc
PUBLIC_HEADER include/openmc.h
LINKER_LANGUAGE Fortran)
target_include_directories(libopenmc
@ -496,8 +494,8 @@ install(TARGETS openmc libopenmc
RUNTIME DESTINATION bin
LIBRARY 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/openmc" RENAME copyright)
install(DIRECTORY include/ DESTINATION include)

View file

@ -1,6 +1,7 @@
include CMakeLists.txt
include LICENSE
include schemas.xml
include pyproject.toml
include openmc/data/reconstruct.pyx
include docs/source/_templates/layout.html
include docs/sphinxext/LICENSE

View file

@ -5,7 +5,7 @@ Windowed Multipole Library Format
=================================
**/version** (*char[]*)
The format version of the file. The current version is "v0.2"
The format version of the file. The current version is "v1.0"
**/nuclide/**
- **broaden_poly** (*int[]*)
@ -23,55 +23,25 @@ Windowed Multipole Library Format
\text{data}[:,i] = [\text{pole},~\text{residue}_1,~\text{residue}_2,
~\ldots]
The residues are in the order: total, competitive if present,
absorption, fission. Complex numbers are stored by forming a type with
":math:`r`" and ":math:`i`" identifiers, similar to how `h5py`_ does it.
- **end_E** (*double*)
The residues are in the order: scattering, absorption, fission. Complex
numbers are stored by forming a type with ":math:`r`" and ":math:`i`"
identifiers, similar to how `h5py`_ does it.
- **E_max** (*double*)
Highest energy the windowed multipole part of the library is valid for.
- **formalism** (*int*)
The formalism of the underlying data. Uses the `ENDF-6`_ format
formalism numbers.
.. table:: Table of supported formalisms.
+-------------+------------------+
| Formalism | Formalism number |
+=============+==================+
| MLBW | 2 |
+-------------+------------------+
| Reich-Moore | 3 |
+-------------+------------------+
- **l_value** (*int[]*)
The index for a corresponding pole. Equivalent to the :math:`l` quantum
number of the resonance the pole comes from :math:`+1`.
- **pseudo_K0RS** (*double[]*)
:math:`l` dependent value of
.. math::
\sqrt{\frac{2 m_n}{\hbar}}\frac{AWR}{AWR + 1} r_{s,l}
Where :math:`m_n` is mass of neutron, :math:`AWR` is the atomic weight
ratio of the target to the neutron, and :math:`r_{s,l}` is the
scattering radius for a given :math:`l`.
- **E_min** (*double*)
Lowest energy the windowed multipole part of the library is valid for.
- **spacing** (*double*)
.. math::
\frac{\sqrt{E_{max}}- \sqrt{E_{min}}}{n_w}
\frac{\sqrt{E_{max}} - \sqrt{E_{min}}}{n_w}
Where :math:`E_{max}` is the maximum energy the windows go up to. This
is not equivalent to the maximum energy for which the windowed multipole
data is valid for. It is slightly higher to ensure an integer number of
windows. :math:`E_{min}` is the minimum energy and equivalent to
``start_E``, and :math:`n_w` is the number of windows, given by
``windows``.
Where :math:`E_{max}` is the maximum energy the windows go up to.
:math:`E_{min}` is the minimum energy, and :math:`n_w` is the number of
windows, given by ``windows``.
- **sqrtAWR** (*double*)
Square root of the atomic weight ratio.
- **start_E** (*double*)
Lowest energy the windowed multipole part of the library is valid for.
- **w_start** (*int[]*)
The pole to start from for each window.
- **w_end** (*int[]*)
The pole to end at for each window.
- **windows** (*int[][]*)
The poles to start from and end at for each window. windows[i, 0] and
windows[i, 1] are, respectively, the indexes (1-based) of the first and
last pole in window i.
.. _h5py: http://docs.h5py.org/en/latest/
.. _ENDF-6: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf

View file

@ -50,6 +50,12 @@ Each ``material`` element can have the following attributes or sub-elements:
*Default*: ""
:depletable:
Boolean value indicating whether the material is depletable.
:volume:
Volume of the material in cm^3.
:temperature:
An element with no attributes which is used to set the default temperature
of the material in Kelvin.

View file

@ -90,7 +90,7 @@ Assuming free-gas thermal motion, cross sections in the multipole form can be
analytically Doppler broadened to give the form:
.. math::
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[i r_j
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[r_j
\sqrt{\pi} W_i(z) - \frac{r_j}{\sqrt{\pi}} C \left(\frac{p_j}{\sqrt{\xi}},
\frac{u}{2 \sqrt{\xi}}\right)\right]
.. math::
@ -141,7 +141,7 @@ scattering does not occur in the resolved resonance region. This is usually,
but not always the case. Future library versions may eliminate this issue.
The data format used by OpenMC to represent windowed multipole data is specified
in :ref:`io_data_wmp`.
in :ref:`io_data_wmp` with a publicly available `WMP library`_.
.. _temperature_treatment:
@ -270,6 +270,7 @@ or even isotropic scattering.
https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
.. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381
.. _Josey: http://dx.doi.org/10.1016/j.jcp.2015.08.013
.. _WMP Library: https://github.com/mit-crpg/WMP_Library
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi
.. _NJOY: http://t2.lanl.gov/codes.shtml

View file

@ -178,6 +178,23 @@ Post-processing
openmc.StatePoint
openmc.Summary
The following classes and functions are used for functional expansion reconstruction.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.ZernikeRadial
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.legendre_from_expcoef
Various classes may be created when performing tally slicing and/or arithmetic:
.. autosummary::

View file

@ -1466,7 +1466,7 @@
},
"outputs": [],
"source": [
"url = 'https://anl.box.com/shared/static/ulhcoohm12gduwdalknmf8dpnepzkxj0.h5'\n",
"url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/092238.h5'\n",
"filename, headers = urllib.request.urlretrieve(url, '092238.h5')"
]
},
@ -1485,7 +1485,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"The `WindowedMultipole` object can be called with energy and temperature values. Calling the object gives a tuple of 3 cross sections: total, radiative capture, and fission."
"The `WindowedMultipole` object can be called with energy and temperature values. Calling the object gives a tuple of 3 cross sections: elastic scattering, radiative capture, and fission."
]
},
{
@ -1498,9 +1498,7 @@
{
"data": {
"text/plain": [
"(array(9.638243132516015),\n",
" array(0.5053244245010787),\n",
" array(2.931753364280356e-06))"
"(array(9.13284265), array(0.50530278), array(2.9316765e-06))"
]
},
"execution_count": 43,

View file

@ -1,5 +1,5 @@
#ifndef OPENMC_H
#define OPENMC_H
#ifndef OPENMC_CAPI_H
#define OPENMC_CAPI_H
#include <stdint.h>
#include <stdbool.h>
@ -37,7 +37,7 @@ extern "C" {
int openmc_filter_set_id(int32_t index, int32_t id);
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_find_cell(double* xyz, int32_t* index, 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);
@ -47,6 +47,7 @@ extern "C" {
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_get_tally_next_id(int32_t* id);
int openmc_hard_reset();
int openmc_init(int argc, char* argv[], const void* intracomm);
int openmc_init_f(const int* intracomm);
@ -56,9 +57,11 @@ extern "C" {
int openmc_material_add_nuclide(int32_t index, const char name[], double density);
int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n);
int openmc_material_get_id(int32_t index, int32_t* id);
int openmc_material_get_volume(int32_t index, double* volume);
int openmc_material_set_density(int32_t index, double density);
int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density);
int openmc_material_set_id(int32_t index, int32_t id);
int openmc_material_set_volume(int32_t index, double volume);
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);
@ -92,15 +95,19 @@ extern "C" {
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
int openmc_statepoint_write(const char filename[]);
int openmc_tally_allocate(int32_t index, const char* type);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
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_get_type(int32_t index, int32_t* type);
int openmc_tally_reset(int32_t index);
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_estimator(int32_t index, const char* estimator);
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);
@ -159,14 +166,14 @@ extern "C" {
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};
const int RUN_MODE_FIXEDSOURCE = 1;
const int RUN_MODE_EIGENVALUE = 2;
const int RUN_MODE_PLOTTING = 3;
const int RUN_MODE_PARTICLE = 4;
const int RUN_MODE_VOLUME = 5;
#ifdef __cplusplus
}
#endif
#endif // OPENMC_H
#endif // OPENMC_CAPI_H

View file

@ -11,7 +11,7 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "position.h"
#include "openmc/position.h"
namespace openmc {

View file

@ -130,17 +130,6 @@ constexpr int ANGLE_HISTOGRAM {5};
constexpr int TEMPERATURE_NEAREST {1};
constexpr int TEMPERATURE_INTERPOLATION {2};
// Secondary energy mode for S(a,b) inelastic scattering
// TODO: Convert to enum
constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins
constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins
constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation
// Elastic mode for S(a,b) elastic scattering
// TODO: Convert to enum
constexpr int SAB_ELASTIC_DISCRETE {3}; // Sample from discrete cosines
constexpr int SAB_ELASTIC_EXACT {4}; // Exact treatment for coherent elastic
// Reaction types
// TODO: Convert to enum
constexpr int TOTAL_XS {1};
@ -260,6 +249,8 @@ constexpr int N_AC {849};
constexpr int N_2N0 {875};
constexpr int N_2NC {891};
constexpr std::array<int, 6> DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N};
// Fission neutron emission (nu) type
constexpr int NU_NONE {0}; // No nu values (non-fissionable)
constexpr int NU_POLYNOMIAL {1}; // Nu values given by polynomial

View file

@ -10,7 +10,7 @@
#include "pugixml.hpp"
#include "constants.h"
#include "openmc/constants.h"
namespace openmc {
@ -107,6 +107,10 @@ public:
//! Sample a value from the distribution
//! \return Sampled value
double sample() const;
// x property
std::vector<double>& x() { return x_; }
const std::vector<double>& x() const { return x_; }
private:
std::vector<double> x_; //!< tabulated independent variable
std::vector<double> p_; //!< tabulated probability density

View file

@ -6,9 +6,10 @@
#include <vector> // for vector
#include "distribution.h"
#include "hdf5.h"
#include "openmc/distribution.h"
namespace openmc {
//==============================================================================

View file

@ -9,8 +9,8 @@
#include "xtensor/xtensor.hpp"
#include "hdf5.h"
#include "constants.h"
#include "endf.h"
#include "openmc/constants.h"
#include "openmc/endf.h"
namespace openmc {

View file

@ -3,8 +3,8 @@
#include <memory>
#include "distribution.h"
#include "position.h"
#include "openmc/distribution.h"
#include "openmc/position.h"
namespace openmc {

View file

@ -3,8 +3,8 @@
#include "pugixml.hpp"
#include "distribution.h"
#include "position.h"
#include "openmc/distribution.h"
#include "openmc/position.h"
namespace openmc {

View file

@ -6,9 +6,10 @@
#include <vector>
#include "constants.h"
#include "hdf5.h"
#include "openmc/constants.h"
namespace openmc {
//! Convert integer representing interpolation law to enum
@ -73,6 +74,19 @@ private:
std::vector<double> y_; //!< values of ordinate
};
//==============================================================================
//! Coherent elastic scattering data from a crystalline material
//==============================================================================
class CoherentElasticXS : public Function1D {
explicit CoherentElasticXS(hid_t dset);
double operator()(double E) const;
private:
std::vector<double> bragg_edges_; //!< Bragg edges in [eV]
std::vector<double> factors_; //!< Partial sums of structure factors [eV-b]
};
} // namespace openmc
#endif // OPENMC_ENDF_H

View file

@ -1,10 +1,11 @@
#ifndef ERROR_H
#define ERROR_H
#ifndef OPENMC_ERROR_H
#define OPENMC_ERROR_H
#include <cstring>
#include <string>
#include <sstream>
#include "openmc/capi.h"
namespace openmc {
@ -14,20 +15,38 @@ extern "C" void warning_from_c(const char* message, int message_len);
extern "C" void write_message_from_c(const char* message, int message_len,
int level);
inline
void fatal_error(const char *message)
inline void
set_errmsg(const char* message)
{
fatal_error_from_c(message, strlen(message));
std::strcpy(openmc_err_msg, message);
}
inline void
set_errmsg(const std::string& message)
{
std::strcpy(openmc_err_msg, message.c_str());
}
inline void
set_errmsg(const std::stringstream& message)
{
std::strcpy(openmc_err_msg, message.str().c_str());
}
inline
void fatal_error(const std::string &message)
void fatal_error(const char* message)
{
fatal_error_from_c(message, std::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)
void fatal_error(const std::stringstream& message)
{
fatal_error(message.str());
}
@ -47,7 +66,7 @@ void warning(const std::stringstream& message)
inline
void write_message(const char* message, int level)
{
write_message_from_c(message, strlen(message), level);
write_message_from_c(message, std::strlen(message), level);
}
inline
@ -63,4 +82,4 @@ void write_message(const std::stringstream& message, int level)
}
} // namespace openmc
#endif // ERROR_H
#endif // OPENMC_ERROR_H

View file

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

View file

@ -1,8 +1,8 @@
//! \file geometry_aux.h
//! Auxilary functions for geometry initialization and general data handling.
#ifndef GEOMETRY_AUX_H
#define GEOMETRY_AUX_H
#ifndef OPENMC_GEOMETRY_AUX_H
#define OPENMC_GEOMETRY_AUX_H
#include <cstdint>
@ -112,4 +112,4 @@ extern "C" int maximum_levels(int32_t univ);
extern "C" void free_memory_geometry_c();
} // namespace openmc
#endif // GEOMETRY_AUX_H
#endif // OPENMC_GEOMETRY_AUX_H

View file

@ -13,7 +13,7 @@
#include "xtensor/xadapt.hpp"
#include "xtensor/xarray.hpp"
#include "position.h"
#include "openmc/position.h"
namespace openmc {
@ -80,6 +80,7 @@ read_nd_vector(hid_t obj_id, const char* name,
bool must_have = false);
std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
std::vector<std::string> dataset_names(hid_t group_id);
void ensure_exists(hid_t group_id, const char* name);
std::vector<std::string> group_names(hid_t group_id);
std::vector<hsize_t> object_shape(hid_t obj_id);
@ -215,12 +216,31 @@ read_attribute(hid_t obj_id, const char* name, std::string& str)
str = std::string{buffer, n};
}
// overload for std::vector<std::string>
inline void
read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
{
auto dims = attribute_shape(obj_id, name);
auto m = dims[0];
// Allocate a C char array to get strings
auto n = attribute_typesize(obj_id, name);
char buffer[m][n+1];
// Read char data in attribute
read_attr_string(obj_id, name, n, buffer[0]);
for (int i = 0; i < m; ++i) {
vec.emplace_back(&buffer[i][0]);
}
}
//==============================================================================
// Templates/overloads for read_dataset
//==============================================================================
template<typename T>
void read_dataset(hid_t obj_id, const char* name, T buffer, bool indep=false)
void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
{
read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
}

View file

@ -1,5 +1,5 @@
#ifndef INITIALIZE_H
#define INITIALIZE_H
#ifndef OPENMC_INITIALIZE_H
#define OPENMC_INITIALIZE_H
#ifdef OPENMC_MPI
#include "mpi.h"
@ -17,4 +17,4 @@ void initialize_mpi(MPI_Comm intracomm);
}
#endif // INITIALIZE_H
#endif // OPENMC_INITIALIZE_H

View file

@ -10,8 +10,8 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "constants.h"
#include "position.h"
#include "openmc/constants.h"
#include "openmc/position.h"
namespace openmc {

View file

@ -25,11 +25,12 @@ class Material
{
public:
int32_t id; //!< Unique ID
double volume_ {-1.0}; //!< Volume in [cm^3]
//! \brief Default temperature for cells containing this material.
//!
//! A negative value indicates no default temperature was specified.
double temperature {-1};
double temperature_ {-1};
Material() {};

View file

@ -7,9 +7,9 @@
#include <cmath>
#include <cstdlib>
#include "constants.h"
#include "position.h"
#include "random_lcg.h"
#include "openmc/constants.h"
#include "openmc/position.h"
#include "openmc/random_lcg.h"
namespace openmc {
@ -93,13 +93,18 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
//==============================================================================
//! Calculate only the even order components of n-th order modified Zernike
//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho)
//! location on the unit disk.
//! Calculate only the even radial components of n-th order modified Zernike
//! polynomial moment with azimuthal dependency m = 0 for a given angle
//! (rho, theta) location on the unit disk.
//!
//! Since m = 0, n could only be even orders. Z_q0 = R_q0
//!
//! See calc_zn_c for methodology.
//! This procedure uses the modified Kintner's method for calculating Zernike
//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
//! R. (2003). A comparative analysis of algorithms for fast computation of
//! Zernike moments. Pattern Recognition, 36(3), 731-742.
//! The normalization of the polynomials is such that the integral of Z_pq^2
//! over the unit disk is exactly pi.
//!
//! @param n The maximum order requested
//! @param rho The radial parameter to specify location on the unit disk

View file

@ -1,5 +1,5 @@
#ifndef MESSAGE_PASSING_H
#define MESSAGE_PASSING_H
#ifndef OPENMC_MESSAGE_PASSING_H
#define OPENMC_MESSAGE_PASSING_H
#ifdef OPENMC_MPI
#include "mpi.h"
@ -19,4 +19,4 @@ namespace mpi {
} // namespace mpi
} // namespace openmc
#endif // MESSAGE_PASSING_H
#endif // OPENMC_MESSAGE_PASSING_H

View file

@ -1,15 +1,15 @@
//! \file mgxs.h
//! A collection of classes for Multi-Group Cross Section data
#ifndef MGXS_H
#define MGXS_H
#ifndef OPENMC_MGXS_H
#define OPENMC_MGXS_H
#include <string>
#include <vector>
#include "constants.h"
#include "hdf5_interface.h"
#include "xsdata.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xsdata.h"
namespace openmc {
@ -202,4 +202,4 @@ class Mgxs {
};
} // namespace openmc
#endif // MGXS_H
#endif // OPENMC_MGXS_H

View file

@ -1,8 +1,8 @@
//! \file mgxs_interface.h
//! A collection of C interfaces to the C++ Mgxs class
#ifndef MGXS_INTERFACE_H
#define MGXS_INTERFACE_H
#ifndef OPENMC_MGXS_INTERFACE_H
#define OPENMC_MGXS_INTERFACE_H
#include "hdf5_interface.h"
#include "mgxs.h"
@ -76,4 +76,4 @@ extern "C" double
get_awr_c(int index);
} // namespace openmc
#endif // MGXS_INTERFACE_H
#endif // OPENMC_MGXS_INTERFACE_H

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

@ -0,0 +1,67 @@
#ifndef OPENMC_NUCLIDE_H
#define OPENMC_NUCLIDE_H
#include "openmc/constants.h"
namespace openmc {
//===============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy
//===============================================================================
struct NuclideMicroXS {
// Microscopic cross sections in barns
double total; //!< total cross section
double absorption; //!< absorption (disappearance)
double fission; //!< fission
double nu_fission; //!< neutron production from fission
double elastic; //!< If sab_frac is not 1 or 0, then this value is
//!< averaged over bound and non-bound nuclei
double thermal; //!< Bound thermal elastic & inelastic scattering
double thermal_elastic; //!< Bound thermal elastic scattering
double photon_prod; //!< microscopic photon production xs
// Cross sections for depletion reactions (note that these are not stored in
// macroscopic cache)
double reaction[DEPLETION_RX.size()];
// Indicies and factors needed to compute cross sections from the data tables
int index_grid; //!< Index on nuclide energy grid
int index_temp; //!< Temperature index for nuclide
double interp_factor; //!< Interpolation factor on nuc. energy grid
int index_sab {-1}; //!< Index in sab_tables
int index_temp_sab; //!< Temperature index for sab_tables
double sab_frac; //!< Fraction of atoms affected by S(a,b)
bool use_ptable; //!< In URR range with probability tables?
// Energy and temperature last used to evaluate these cross sections. If
// these values have changed, then the cross sections must be re-evaluated.
double last_E {0.0}; //!< Last evaluated energy
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
//!< * temperature (eV))
};
//===============================================================================
// MATERIALMACROXS contains cached macroscopic cross sections for the material a
// particle is traveling through
//===============================================================================
struct MaterialMacroXS {
double total; //!< macroscopic total xs
double absorption; //!< macroscopic absorption xs
double fission; //!< macroscopic fission xs
double nu_fission; //!< macroscopic production xs
double photon_prod; //!< macroscopic photon production xs
// Photon cross sections
double coherent; //!< macroscopic coherent xs
double incoherent; //!< macroscopic incoherent xs
double photoelectric; //!< macroscopic photoelectric xs
double pair_production; //!< macroscopic pair production xs
};
} // namespace openmc
#endif // OPENMC_NUCLIDE_H

View file

@ -9,7 +9,7 @@
#include <sstream>
#include <string>
#include "openmc.h"
#include "openmc/capi.h"
namespace openmc {

View file

@ -1,5 +1,5 @@
#ifndef PLOT_H
#define PLOT_H
#ifndef OPENMC_PLOT_H
#define OPENMC_PLOT_H
#include "hdf5.h"
@ -12,4 +12,4 @@ extern "C" void voxel_write_slice(int x, hid_t dspace, hid_t dset,
extern "C" void voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace);
} // namespace openmc
#endif // PLOT_H
#endif // OPENMC_PLOT_H

View file

@ -1,5 +1,5 @@
#ifndef RANDOM_LCG_H
#define RANDOM_LCG_H
#ifndef OPENMC_RANDOM_LCG_H
#define OPENMC_RANDOM_LCG_H
#include <cstdint>
@ -93,4 +93,4 @@ extern "C" int64_t openmc_get_seed();
extern "C" void openmc_set_seed(int64_t new_seed);
} // namespace openmc
#endif // RANDOM_LCG_H
#endif // OPENMC_RANDOM_LCG_H

View file

@ -7,7 +7,8 @@
#include <vector>
#include "hdf5.h"
#include "reaction_product.h"
#include "openmc/reaction_product.h"
namespace openmc {

View file

@ -8,9 +8,10 @@
#include <vector> // for vector
#include "hdf5.h"
#include "angle_energy.h"
#include "endf.h"
#include "particle.h"
#include "openmc/angle_energy.h"
#include "openmc/endf.h"
#include "openmc/particle.h"
namespace openmc {

View file

@ -1,11 +1,13 @@
//! \file scattdata.h
//! A collection of multi-group scattering data classes
#ifndef SCATTDATA_H
#define SCATTDATA_H
#ifndef OPENMC_SCATTDATA_H
#define OPENMC_SCATTDATA_H
#include <vector>
#include "openmc/constants.h"
namespace openmc {
// forward declarations so we can name our friend functions
@ -257,4 +259,4 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu);
} // namespace openmc
#endif // SCATTDATA_H
#endif // OPENMC_SCATTDATA_H

View file

@ -8,9 +8,10 @@
#include "hdf5.h"
#include "xtensor/xtensor.hpp"
#include "angle_energy.h"
#include "endf.h"
#include "distribution.h"
#include "openmc/angle_energy.h"
#include "openmc/endf.h"
#include "openmc/distribution.h"
namespace openmc {
@ -21,14 +22,6 @@ namespace openmc {
class CorrelatedAngleEnergy : public AngleEnergy {
public:
explicit CorrelatedAngleEnergy(hid_t group);
//! Sample distribution for an angle and energy
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
void sample(double E_in, double& E_out, double& mu) const;
private:
//! Outgoing energy/angle at a single incoming energy
struct CorrTable {
int n_discrete; //!< Number of discrete lines
@ -39,6 +32,22 @@ private:
std::vector<UPtrDist> angle; //!< Angle distribution
};
explicit CorrelatedAngleEnergy(hid_t group);
//! Sample distribution for an angle and energy
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
void sample(double E_in, double& E_out, double& mu) const;
// energy property
std::vector<double>& energy() { return energy_; }
const std::vector<double>& energy() const { return energy_; }
// distribution property
std::vector<CorrTable>& distribution() { return distribution_; }
const std::vector<CorrTable>& distribution() const { return distribution_; }
private:
int n_region_; //!< Number of interpolation regions
std::vector<int> breakpoints_; //!< Breakpoints between regions
std::vector<Interpolation> interpolation_; //!< Interpolation laws

View file

@ -8,9 +8,10 @@
#include "hdf5.h"
#include "xtensor/xtensor.hpp"
#include "angle_energy.h"
#include "constants.h"
#include "endf.h"
#include "openmc/angle_energy.h"
#include "openmc/constants.h"
#include "openmc/endf.h"
namespace openmc {

View file

@ -6,7 +6,7 @@
#include "hdf5.h"
#include "angle_energy.h"
#include "openmc/angle_energy.h"
namespace openmc {

View file

@ -8,9 +8,10 @@
#include <vector>
#include "hdf5.h"
#include "angle_energy.h"
#include "distribution_angle.h"
#include "distribution_energy.h"
#include "openmc/angle_energy.h"
#include "openmc/distribution_angle.h"
#include "openmc/distribution_energy.h"
namespace openmc {

View file

@ -4,6 +4,7 @@
//! \file settings.h
//! \brief Settings for OpenMC
#include <array>
#include <string>
#include "pugixml.hpp"
@ -35,6 +36,12 @@ extern std::string path_multipole;
extern std::string path_output;
extern std::string path_source;
extern int temperature_method;
extern bool temperature_multipole;
extern double temperature_tolerance;
extern double temperature_default;
extern std::array<double, 2> temperature_range;
//==============================================================================
//! Read settings from XML file
//! \param[in] root XML node for <settings>

View file

@ -1,5 +1,5 @@
#ifndef SIMULATION_H
#define SIMULATION_H
#ifndef OPENMC_SIMULATION_H
#define OPENMC_SIMULATION_H
#include <cstdint>
@ -10,4 +10,4 @@ extern "C" int openmc_n_lost_particles;
#pragma omp threadprivate(openmc_current_work)
#endif // SIMULATION_H
#endif // OPENMC_SIMULATION_H

View file

@ -1,10 +1,11 @@
#ifndef STATE_POINT_H
#define STATE_POINT_H
#ifndef OPENMC_STATE_POINT_H
#define OPENMC_STATE_POINT_H
#include <cstdint>
#include "hdf5.h"
#include "openmc.h"
#include "openmc/capi.h"
namespace openmc {
@ -14,4 +15,4 @@ extern "C" void read_source_bank(hid_t group_id, int64_t* work_index,
Bank* source_bank);
} // namespace openmc
#endif // STATE_POINT_H
#endif // OPENMC_STATE_POINT_H

View file

@ -1,8 +1,8 @@
//! \file string_functions.h
//! A collection of helper routines for C-strings and STL strings
#ifndef STRING_FUNCTIONS_H
#define STRING_FUNCTIONS_H
#ifndef OPENMC_STRING_FUNCTIONS_H
#define OPENMC_STRING_FUNCTIONS_H
#include <string>

View file

@ -9,8 +9,8 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "constants.h"
#include "position.h"
#include "openmc/constants.h"
#include "openmc/position.h"
namespace openmc {

147
include/openmc/thermal.h Normal file
View file

@ -0,0 +1,147 @@
#ifndef OPENMC_THERMAL_H
#define OPENMC_THERMAL_H
#include <cstddef>
#include <string>
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/hdf5_interface.h"
#include "openmc/nuclide.h"
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
// Secondary energy mode for S(a,b) inelastic scattering
// TODO: Convert to enum
constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins
constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins
constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation
// Elastic mode for S(a,b) elastic scattering
// TODO: Convert to enum
constexpr int SAB_ELASTIC_INCOHERENT {3}; // Incoherent elastic scattering
constexpr int SAB_ELASTIC_COHERENT {4}; // Coherent elastic scattering (Bragg edges)
//==============================================================================
//! Secondary angle-energy data for thermal neutron scattering at a single
//! temperature
//==============================================================================
class ThermalData {
public:
ThermalData(hid_t group, int secondary_mode);
// Sample an outgoing energy and angle
void sample(const NuclideMicroXS* micro_xs, double E_in,
double* E_out, double* mu);
private:
//! Secondary energy/angle distributions for inelastic thermal scattering
//! collisions which utilize a continuous secondary energy representation.
struct DistEnergySab {
std::size_t n_e_out; //!< Number of outgoing energies
xt::xtensor<double, 1> e_out; //!< Outgoing energies
xt::xtensor<double, 1> e_out_pdf; //!< Probability density function
xt::xtensor<double, 1> e_out_cdf; //!< Cumulative distribution function
xt::xtensor<double, 2> mu; //!< Equiprobable angles at each outgoing energy
};
//! Upper threshold for incoherent inelastic scattering (usually ~4 eV)
double threshold_inelastic_;
//! Upper threshold for coherent/incoherent elastic scattering
double threshold_elastic_ {0.0};
// Inelastic scattering data
int inelastic_mode_; //!< distribution type (equal/skewed/continuous)
std::size_t n_inelastic_e_in_; //!< number of incoming E for inelastic
std::size_t n_inelastic_e_out_; //!< number of outgoing E for inelastic
std::size_t n_inelastic_mu_; //!< number of outgoing angles for inelastic
std::vector<double> inelastic_e_in_; //!< incoming E grid for inelastic
std::vector<double> inelastic_sigma_; //!< inelastic scattering cross section
// The following are used only for equal/skewed distributions
xt::xtensor<double, 2> inelastic_e_out_;
xt::xtensor<double, 3> inelastic_mu_;
// The following is used only for continuous S(a,b) distributions. The
// different implementation is necessary because the continuous representation
// has a variable number of outgoing energy points for each incoming energy
std::vector<DistEnergySab> inelastic_data_; //!< Secondary angle-energy at
//!< each incoming energy
// Elastic scattering data
int elastic_mode_; //!< type of elastic (incoherent/coherent)
std::size_t n_elastic_e_in_; //!< number of incoming E for elastic
std::size_t n_elastic_mu_; //!< number of outgoing angles for elastic
std::vector<double> elastic_e_in_; //!< incoming E grid for elastic
std::vector<double> elastic_P_; //!< elastic scattering cross section
xt::xtensor<double, 2> elastic_mu_; //!< equi-probable angles at each incoming E
// ThermalScattering needs access to private data members
friend class ThermalScattering;
};
//==============================================================================
//! Data for thermal neutron scattering, typically off light isotopes in
//! moderating materials such as water, graphite, BeO, etc.
//==============================================================================
class ThermalScattering {
public:
ThermalScattering(hid_t group, const std::vector<double>& temperature, int method,
double tolerance, const double* minmax);
//! Determine inelastic/elastic cross section at given energy
//!
//! \param[in] E incoming energy in [eV]
//! \param[in] sqrtkT square-root of temperature multipled by Boltzmann's constant
//! \param[out] i_temp corresponding temperature index
//! \param[out] elastic Thermal elastic scattering cross section
//! \param[out] inelastic Thermal inelastic scattering cross section
void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
double* inelastic) const;
//! Determine whether table applies to a particular nuclide
//!
//! \param[in] name Name of the nuclide, e.g., "H1"
//! \return Whether table applies to the nuclide
bool has_nuclide(const char* name) const;
// Sample an outgoing energy and angle
void sample(const NuclideMicroXS* micro_xs, double E_in,
double* E_out, double* mu);
double threshold() const { return data_[0].threshold_inelastic_; }
std::string name_; //!< name of table, e.g. "c_H_in_H2O"
double awr_; //!< weight of nucleus in neutron masses
std::vector<double> kTs_; //!< temperatures in [eV] (k*T)
std::vector<std::string> nuclides_; //!< Valid nuclides
//! cross sections and distributions at each temperature
std::vector<ThermalData> data_;
};
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature,
int n, int method, double tolerance, const double* minmax);
void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT,
int* i_temp, double* elastic, double* inelastic);
void sab_free(ThermalScattering* data);
bool sab_has_nuclide(ThermalScattering* data, const char* name);
void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs,
double E_in, double* E_out, double* mu);
double sab_threshold(ThermalScattering* data);
}
} // namespace openmc
#endif // OPENMC_THERMAL_H

View file

@ -1,5 +1,5 @@
#ifndef XML_INTERFACE_H
#define XML_INTERFACE_H
#ifndef OPENMC_XML_INTERFACE_H
#define OPENMC_XML_INTERFACE_H
#include <sstream> // for stringstream
#include <string>
@ -16,9 +16,11 @@ check_for_node(pugi::xml_node node, const char *name)
return node.attribute(name) || node.child(name);
}
std::string get_node_value(pugi::xml_node node, const char *name,
std::string get_node_value(pugi::xml_node node, const char* name,
bool lowercase=false, bool strip=false);
bool get_node_value_bool(pugi::xml_node node, const char* name);
template <typename T>
std::vector<T> get_node_array(pugi::xml_node node, const char* name,
bool lowercase=false)
@ -37,4 +39,4 @@ std::vector<T> get_node_array(pugi::xml_node node, const char* name,
}
} // namespace openmc
#endif // XML_INTERFACE_H
#endif // OPENMC_XML_INTERFACE_H

View file

@ -1,14 +1,14 @@
//! \file xsdata.h
//! A collection of classes for containing the Multi-Group Cross Section data
#ifndef XSDATA_H
#define XSDATA_H
#ifndef OPENMC_XSDATA_H
#define OPENMC_XSDATA_H
#include <memory>
#include <vector>
#include "hdf5_interface.h"
#include "scattdata.h"
#include "openmc/hdf5_interface.h"
#include "openmc/scattdata.h"
namespace openmc {
@ -115,4 +115,4 @@ class XsData {
} //namespace openmc
#endif // XSDATA_H
#endif // OPENMC_XSDATA_H

View file

@ -28,6 +28,7 @@ from openmc.particle_restart import *
from openmc.mixin import *
from openmc.plotter import *
from openmc.search import *
from openmc.polynomial import *
from . import examples
# Import a few convencience functions that used to be here

51
openmc/_xml.py Normal file
View file

@ -0,0 +1,51 @@
def clean_indentation(element, level=0, spaces_per_level=2):
"""
copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is
printed in a nice way
"""
i = "\n" + level*spaces_per_level*" "
if len(element):
if not element.text or not element.text.strip():
element.text = i + spaces_per_level*" "
if not element.tail or not element.tail.strip():
element.tail = i
for sub_element in element:
clean_indentation(sub_element, level+1, spaces_per_level)
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):
element.tail = i
def get_text(elem, name, default=None):
"""Retrieve text of an attribute or subelement.
Parameters
----------
elem : xml.etree.ElementTree.Element
Element from which to search
name : str
Name of attribute/subelement
default : object
A defult value to return if matching attribute/subelement exists
Returns
-------
str
Text of attribute or subelement
"""
if name in elem.attrib:
return elem.get(name, default)
else:
child = elem.find(name)
return child.text if child is not None else default

View file

@ -24,10 +24,10 @@ _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_find_cell.argtypes = [POINTER(c_double*3), POINTER(c_int32),
POINTER(c_int32)]
_dll.openmc_find_cell.restype = c_int
_dll.openmc_find_cell.errcheck = _error_handler
_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]
@ -84,10 +84,10 @@ def find_cell(xyz):
indicates which instance it is, i.e., 0 would be the first instance.
"""
uid = c_int32()
index = c_int32()
instance = c_int32()
_dll.openmc_find((c_double*3)(*xyz), 1, uid, instance)
return openmc.capi.cells[uid.value], instance.value
_dll.openmc_find_cell((c_double*3)(*xyz), index, instance)
return openmc.capi.Cell(index=index.value), instance.value
def find_material(xyz):
@ -104,11 +104,15 @@ def find_material(xyz):
Material containing the point, or None is no material is found
"""
uid = c_int32()
index = c_int32()
instance = c_int32()
_dll.openmc_find((c_double*3)(*xyz), 2, uid, instance)
return openmc.capi.materials[uid.value] if uid != 0 else None
_dll.openmc_find_cell((c_double*3)(*xyz), index, instance)
mats = openmc.capi.Cell(index=index.value).fill
if isinstance(mats, openmc.capi.Material):
return mats
else:
return mats[instance]
def hard_reset():
"""Reset tallies, timers, and pseudo-random number generator state."""

View file

@ -5,7 +5,7 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from openmc.exceptions import AllocationError, InvalidIDError
from openmc.exceptions import AllocationError, InvalidIDError, OpenMCError
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler
@ -32,6 +32,9 @@ _dll.openmc_material_get_densities.argtypes = [
POINTER(c_int)]
_dll.openmc_material_get_densities.restype = c_int
_dll.openmc_material_get_densities.errcheck = _error_handler
_dll.openmc_material_get_volume.argtypes = [c_int32, POINTER(c_double)]
_dll.openmc_material_get_volume.restype = c_int
_dll.openmc_material_get_volume.errcheck = _error_handler
_dll.openmc_material_set_density.argtypes = [c_int32, c_double]
_dll.openmc_material_set_density.restype = c_int
_dll.openmc_material_set_density.errcheck = _error_handler
@ -42,6 +45,9 @@ _dll.openmc_material_set_densities.errcheck = _error_handler
_dll.openmc_material_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_material_set_id.restype = c_int
_dll.openmc_material_set_id.errcheck = _error_handler
_dll.openmc_material_set_volume.argtypes = [c_int32, c_double]
_dll.openmc_material_set_volume.restype = c_int
_dll.openmc_material_set_volume.errcheck = _error_handler
class Material(_FortranObjectWithID):
@ -113,6 +119,19 @@ class Material(_FortranObjectWithID):
def id(self, mat_id):
_dll.openmc_material_set_id(self._index, mat_id)
@property
def volume(self):
volume = c_double()
try:
_dll.openmc_material_get_volume(self._index, volume)
except OpenMCError:
return None
return volume.value
@volume.setter
def volume(self, volume):
_dll.openmc_material_set_volume(self._index, volume)
@property
def nuclides(self):
return self._get_densities()[0]

View file

@ -26,9 +26,15 @@ _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_allocate.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_allocate.restype = c_int
_dll.openmc_tally_allocate.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_estimator.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_estimator.restype = c_int
_dll.openmc_tally_get_estimator.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 +53,12 @@ _dll.openmc_tally_get_scores.argtypes = [
c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_tally_get_scores.restype = c_int
_dll.openmc_tally_get_scores.errcheck = _error_handler
_dll.openmc_tally_get_type.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_type.restype = c_int
_dll.openmc_tally_get_type.errcheck = _error_handler
_dll.openmc_tally_reset.argtypes = [c_int32]
_dll.openmc_tally_reset.restype = c_int
_dll.openmc_tally_reset.errcheck = _error_handler
_dll.openmc_tally_results.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)]
_dll.openmc_tally_results.restype = c_int
@ -57,6 +69,9 @@ _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
_dll.openmc_tally_set_estimator.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_set_estimator.restype = c_int
_dll.openmc_tally_set_estimator.errcheck = _error_handler
_dll.openmc_tally_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_tally_set_id.restype = c_int
_dll.openmc_tally_set_id.errcheck = _error_handler
@ -78,6 +93,12 @@ _SCORES = {
-12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-15: 'fission-q-recoverable', -16: 'decay-rate'
}
_ESTIMATORS = {
1: 'analog', 2: 'tracklength', 3: 'collision'
}
_TALLY_TYPES = {
1: 'volume', 2: 'mesh-surface', 3: 'surface'
}
def global_tallies():
@ -140,6 +161,8 @@ class Tally(_FortranObjectWithID):
----------
id : int
ID of the tally
estimator: str
Estimator type of tally (analog, tracklength, collision)
filters : list
List of tally filters
mean : numpy.ndarray
@ -152,6 +175,8 @@ class Tally(_FortranObjectWithID):
Array of tally results
std_dev : numpy.ndarray
An array containing the sample standard deviation for each bin
type : str
Type of tally (volume, mesh_surface, surface)
"""
__instances = WeakValueDictionary()
@ -170,7 +195,7 @@ class Tally(_FortranObjectWithID):
index = c_int32()
_dll.openmc_extend_tallies(1, index, None)
_dll.openmc_tally_set_type(index, b'generic')
_dll.openmc_tally_allocate(index, b'generic')
index = index.value
else:
index = mapping[uid]._index
@ -190,6 +215,26 @@ class Tally(_FortranObjectWithID):
_dll.openmc_tally_get_active(self._index, active)
return active.value
@property
def type(self):
type = c_int32()
_dll.openmc_tally_get_type(self._index, type)
return _TALLY_TYPES[type.value]
@type.setter
def type(self, type):
_dll.openmc_tally_set_type(self._index, type.encode())
@property
def estimator(self):
estimator = c_int32()
_dll.openmc_tally_get_estimator(self._index, estimator)
return _ESTIMATORS[estimator.value]
@estimator.setter
def estimator(self, estimator):
_dll.openmc_tally_set_estimator(self._index, estimator.encode())
@active.setter
def active(self, active):
_dll.openmc_tally_set_active(self._index, active)
@ -302,6 +347,10 @@ class Tally(_FortranObjectWithID):
return std_dev
def reset(self):
"""Reset results and num_realizations of tally"""
_dll.openmc_tally_reset(self._index)
def ci_width(self, alpha=0.05):
"""Confidence interval half-width based on a Student t distribution

View file

@ -13,6 +13,7 @@ import openmc
import openmc.checkvalue as cv
from openmc.surface import Halfspace
from openmc.region import Region, Intersection, Complement
from openmc._xml import get_text
from .mixin import IDManagerMixin
@ -522,3 +523,61 @@ class Cell(IDManagerMixin):
element.set("rotation", ' '.join(map(str, self.rotation)))
return element
@classmethod
def from_xml_element(cls, elem, surfaces, materials, get_universe):
"""Generate cell from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
`<cell>` element
surfaces : dict
Dictionary mapping surface IDs to :class:`openmc.Surface` instances
materials : dict
Dictionary mapping material IDs to :class:`openmc.Material`
instances (defined in :math:`openmc.Geometry.from_xml`)
get_universe : function
Function returning universe (defined in
:meth:`openmc.Geometry.from_xml`)
Returns
-------
Cell
Cell instance
"""
cell_id = int(get_text(elem, 'id'))
name = get_text(elem, 'name')
c = cls(cell_id, name)
# Assign material/distributed materials or fill
mat_text = get_text(elem, 'material')
if mat_text is not None:
mat_ids = mat_text.split()
if len(mat_ids) > 1:
c.fill = [materials[i] for i in mat_ids]
else:
c.fill = materials[mat_ids[0]]
else:
fill_id = int(get_text(elem, 'fill'))
c.fill = get_universe(fill_id)
# Assign region
region = get_text(elem, 'region')
if region is not None:
c.region = Region.from_expression(region, surfaces)
# Check for other attributes
t = get_text(elem, 'temperature')
if t is not None:
c.temperature = float(t)
for key in ('temperature', 'rotation', 'translation'):
value = get_text(elem, key)
if value is not None:
setattr(c, key, [float(x) for x in value.split()])
# Add this cell to appropriate universe
univ_id = int(get_text(elem, 'universe', 0))
get_universe(univ_id).add_cell(c)
return c

View file

@ -1,26 +0,0 @@
def clean_xml_indentation(element, level=0, spaces_per_level=2):
"""
copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
it basically walks your tree and adds spaces and newlines so the tree is
printed in a nice way
"""
i = "\n" + level*spaces_per_level*" "
if len(element):
if not element.text or not element.text.strip():
element.text = i + spaces_per_level*" "
if not element.tail or not element.tail.strip():
element.tail = i
for sub_element in element:
clean_xml_indentation(sub_element, level+1, spaces_per_level)
if not sub_element.tail or not sub_element.tail.strip():
sub_element.tail = i
else:
if level and (not element.tail or not element.tail.strip()):
element.tail = i

View file

@ -15,7 +15,7 @@ from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
@ -512,7 +512,7 @@ class CMFD(object):
self._create_write_matrices_subelement()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._cmfd_file)
clean_indentation(self._cmfd_file)
# Write the XML Tree to the cmfd.xml file
tree = ET.ElementTree(self._cmfd_file)

View file

@ -4,7 +4,7 @@ HDF5_VERSION_MINOR = 0
HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
# Version of WMP nuclear data format
WMP_VERSION = 'v0.2'
WMP_VERSION = 'v1.0'
from .data import *

View file

@ -4,7 +4,7 @@ import xml.etree.ElementTree as ET
import h5py
from openmc.mixin import EqualityMixin
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from openmc.checkvalue import check_type
@ -95,7 +95,7 @@ class DataLibrary(EqualityMixin):
lib_element.set('type', library['type'])
# Clean the indentation to be user-readable
clean_xml_indentation(root)
clean_indentation(root)
# Write XML file
tree = ET.ElementTree(root)

View file

@ -10,26 +10,16 @@ import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
# Formalisms
_FORM_MLBW = 2
_FORM_RM = 3
# Constants that determine which value to access
_MP_EA = 0 # Pole
# Reich-Moore indices
_RM_RT = 1 # Residue total
_RM_RA = 2 # Residue absorption
_RM_RF = 3 # Residue fission
# Multi-level Breit Wigner indices
_MLBW_RT = 1 # Residue total
_MLBW_RX = 2 # Residue competitive
_MLBW_RA = 3 # Residue absorption
_MLBW_RF = 4 # Residue fission
# Residue indices
_MP_RS = 1 # Residue scattering
_MP_RA = 2 # Residue absorption
_MP_RF = 3 # Residue fission
# Polynomial fit indices
_FIT_T = 0 # Total
_FIT_S = 0 # Scattering
_FIT_A = 1 # Absorption
_FIT_F = 2 # Fission
@ -143,98 +133,61 @@ class WindowedMultipole(EqualityMixin):
Parameters
----------
formalism : {'MLBW', 'RM'}
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
for multi-level Breit Wigner or 'RM' for Reich-Moore.
Attributes
----------
num_l : Integral
Number of possible l quantum states for this nuclide.
fit_order : Integral
Order of the windowed curvefit.
fissionable : bool
Whether or not the target nuclide has fission data.
formalism : {'MLBW', 'RM'}
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
for multi-level Breit Wigner or 'RM' for Reich-Moore.
spacing : Real
The width of each window in sqrt(E)-space. For example, the frst window
will end at (sqrt(start_E) + spacing)**2 and the second window at
(sqrt(start_E) + 2*spacing)**2.
will end at (sqrt(E_min) + spacing)**2 and the second window at
(sqrt(E_min) + 2*spacing)**2.
sqrtAWR : Real
Square root of the atomic weight ratio of the target nuclide.
start_E : Real
E_min : Real
Lowest energy in eV the library is valid for.
end_E : Real
E_max : Real
Highest energy in eV the library is valid for.
data : np.ndarray
A 2D array of complex poles and residues. data[i, 0] gives the energy
at which pole i is located. data[i, 1:] gives the residues associated
with the i-th pole. There are 3 residues for Reich-Moore data, one each
for the total, absorption, and fission channels. Multi-level
Breit Wigner data has an additional residue for the competitive channel.
pseudo_k0RS : np.ndarray
A 1D array of Real values. There is one value for each valid l
quantum number. The values are equal to
sqrt(2 m / hbar) * AWR / (AWR + 1) * r
where m is the neutron mass, AWR is the atomic weight ratio, and r
is the l-dependent scattering radius.
l_value : np.ndarray
A 1D array of Integral values equal to the l quantum number for each
pole + 1.
w_start : np.ndarray
A 1D array of Integral values. w_start[i] - 1 is the index of the first
pole in window i.
w_end : np.ndarray
A 1D array of Integral values. w_end[i] - 1 is the index of the last
pole in window i.
with the i-th pole. There are 3 residues, one each for the scattering,
absorption, and fission channels.
windows : np.ndarray
A 2D array of Integral values. windows[i, 0] - 1 is the index of the
first pole in window i. windows[i, 1] - 1 is the index of the last pole
in window i.
broaden_poly : np.ndarray
A 1D array of boolean values indicating whether or not the polynomial
curvefit in that window should be Doppler broadened.
curvefit : np.ndarray
A 3D array of Real curvefit polynomial coefficients. curvefit[i, 0, :]
gives coefficients for the total cross section in window i.
gives coefficients for the scattering cross section in window i.
curvefit[i, 1, :] gives absorption coefficients and curvefit[i, 2, :]
gives fission coefficients. The polynomial terms are increasing powers
of sqrt(E) starting with 1/E e.g:
a/E + b/sqrt(E) + c + d sqrt(E) + ...
"""
def __init__(self, formalism):
self._num_l = None
self.formalism = formalism
def __init__(self):
self.spacing = None
self.sqrtAWR = None
self.start_E = None
self.end_E = None
self.E_min = None
self.E_max = None
self.data = None
self.pseudo_k0RS = None
self.l_value = None
self.w_start = None
self.w_end = None
self.windows = None
self.broaden_poly = None
self.curvefit = None
@property
def num_l(self):
return self._num_l
@property
def fit_order(self):
return self.curvefit.shape[1] - 1
@property
def fissionable(self):
if self.formalism == 'RM':
return self.data.shape[1] == 4
else:
# Assume self.formalism == 'MLBW'
return self.data.shape[1] == 5
@property
def formalism(self):
return self._formalism
return self.data.shape[1] == 4
@property
def spacing(self):
@ -245,32 +198,20 @@ class WindowedMultipole(EqualityMixin):
return self._sqrtAWR
@property
def start_E(self):
return self._start_E
def E_min(self):
return self._E_min
@property
def end_E(self):
return self._end_E
def E_max(self):
return self._E_max
@property
def data(self):
return self._data
@property
def pseudo_k0RS(self):
return self._pseudo_k0RS
@property
def l_value(self):
return self._l_value
@property
def w_start(self):
return self._w_start
@property
def w_end(self):
return self._w_end
def windows(self):
return self._windows
@property
def broaden_poly(self):
@ -280,12 +221,6 @@ class WindowedMultipole(EqualityMixin):
def curvefit(self):
return self._curvefit
@formalism.setter
def formalism(self, formalism):
cv.check_type('formalism', formalism, str)
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
self._formalism = formalism
@spacing.setter
def spacing(self, spacing):
if spacing is not None:
@ -300,19 +235,19 @@ class WindowedMultipole(EqualityMixin):
cv.check_greater_than('sqrtAWR', sqrtAWR, 0.0, equality=False)
self._sqrtAWR = sqrtAWR
@start_E.setter
def start_E(self, start_E):
if start_E is not None:
cv.check_type('start_E', start_E, Real)
cv.check_greater_than('start_E', start_E, 0.0, equality=True)
self._start_E = start_E
@E_min.setter
def E_min(self, E_min):
if E_min is not None:
cv.check_type('E_min', E_min, Real)
cv.check_greater_than('E_min', E_min, 0.0, equality=True)
self._E_min = E_min
@end_E.setter
def end_E(self, end_E):
if end_E is not None:
cv.check_type('end_E', end_E, Real)
cv.check_greater_than('end_E', end_E, 0.0, equality=False)
self._end_E = end_E
@E_max.setter
def E_max(self, E_max):
if E_max is not None:
cv.check_type('E_max', E_max, Real)
cv.check_greater_than('E_max', E_max, 0.0, equality=False)
self._E_max = E_max
@data.setter
def data(self, data):
@ -320,71 +255,25 @@ class WindowedMultipole(EqualityMixin):
cv.check_type('data', data, np.ndarray)
if len(data.shape) != 2:
raise ValueError('Multipole data arrays must be 2D')
if self.formalism == 'RM':
if data.shape[1] not in (3, 4):
raise ValueError('For the Reich-Moore formalism, '
'data.shape[1] must be 3 or 4. One value for the pole.'
' One each for the total and absorption residues. '
'Possibly one more for a fission residue.')
else:
# Assume self.formalism == 'MLBW'
if data.shape[1] not in (4, 5):
raise ValueError('For the Multi-level Breit-Wigner '
'formalism, data.shape[1] must be 4 or 5. One value '
'for the pole. One each for the total, competitive, '
'and absorption residues. Possibly one more for a '
'fission residue.')
if not np.issubdtype(data.dtype, complex):
if data.shape[1] not in (3, 4):
raise ValueError(
'data.shape[1] must be 3 or 4. One value for the pole.'
' One each for the scattering and absorption residues. '
'Possibly one more for a fission residue.')
if not np.issubdtype(data.dtype, np.complexfloating):
raise TypeError('Multipole data arrays must be complex dtype')
self._data = data
@pseudo_k0RS.setter
def pseudo_k0RS(self, pseudo_k0RS):
if pseudo_k0RS is not None:
cv.check_type('pseudo_k0RS', pseudo_k0RS, np.ndarray)
if len(pseudo_k0RS.shape) != 1:
raise ValueError('Multipole pseudo_k0RS arrays must be 1D')
if not np.issubdtype(pseudo_k0RS.dtype, float):
raise TypeError('Multipole data arrays must be float dtype')
self._pseudo_k0RS = pseudo_k0RS
@l_value.setter
def l_value(self, l_value):
if l_value is not None:
cv.check_type('l_value', l_value, np.ndarray)
if len(l_value.shape) != 1:
raise ValueError('Multipole l_value arrays must be 1D')
if not np.issubdtype(l_value.dtype, int):
raise TypeError('Multipole l_value arrays must be integer'
@windows.setter
def windows(self, windows):
if windows is not None:
cv.check_type('windows', windows, np.ndarray)
if len(windows.shape) != 2:
raise ValueError('Multipole windows arrays must be 2D')
if not np.issubdtype(windows.dtype, np.integer):
raise TypeError('Multipole windows arrays must be integer'
' dtype')
self._num_l = len(np.unique(l_value))
else:
self._num_l = None
self._l_value = l_value
@w_start.setter
def w_start(self, w_start):
if w_start is not None:
cv.check_type('w_start', w_start, np.ndarray)
if len(w_start.shape) != 1:
raise ValueError('Multipole w_start arrays must be 1D')
if not np.issubdtype(w_start.dtype, int):
raise TypeError('Multipole w_start arrays must be integer'
' dtype')
self._w_start = w_start
@w_end.setter
def w_end(self, w_end):
if w_end is not None:
cv.check_type('w_end', w_end, np.ndarray)
if len(w_end.shape) != 1:
raise ValueError('Multipole w_end arrays must be 1D')
if not np.issubdtype(w_end.dtype, int):
raise TypeError('Multipole w_end arrays must be integer dtype')
self._w_end = w_end
self._windows = windows
@broaden_poly.setter
def broaden_poly(self, broaden_poly):
@ -392,7 +281,7 @@ class WindowedMultipole(EqualityMixin):
cv.check_type('broaden_poly', broaden_poly, np.ndarray)
if len(broaden_poly.shape) != 1:
raise ValueError('Multipole broaden_poly arrays must be 1D')
if not np.issubdtype(broaden_poly.dtype, bool):
if not np.issubdtype(broaden_poly.dtype, np.bool_):
raise TypeError('Multipole broaden_poly arrays must be boolean'
' dtype')
self._broaden_poly = broaden_poly
@ -403,10 +292,10 @@ class WindowedMultipole(EqualityMixin):
cv.check_type('curvefit', curvefit, np.ndarray)
if len(curvefit.shape) != 3:
raise ValueError('Multipole curvefit arrays must be 3D')
if curvefit.shape[2] not in (2, 3): # sig_t, sig_a (maybe sig_f)
if curvefit.shape[2] not in (2, 3): # sig_s, sig_a (maybe sig_f)
raise ValueError('The third dimension of multipole curvefit'
' arrays must have a length of 2 or 3')
if not np.issubdtype(curvefit.dtype, float):
if not np.issubdtype(curvefit.dtype, np.floating):
raise TypeError('Multipole curvefit arrays must be float dtype')
self._curvefit = curvefit
@ -443,19 +332,14 @@ class WindowedMultipole(EqualityMixin):
'Python API expects version ' + WMP_VERSION)
group = h5file['nuclide']
# Read scalars.
out = cls()
if group['formalism'].value == _FORM_MLBW:
out = cls('MLBW')
elif group['formalism'].value == _FORM_RM:
out = cls('RM')
else:
raise ValueError('Unrecognized/Unsupported R-matrix formalism')
# Read scalars.
out.spacing = group['spacing'].value
out.sqrtAWR = group['sqrtAWR'].value
out.start_E = group['start_E'].value
out.end_E = group['end_E'].value
out.E_min = group['E_min'].value
out.E_max = group['E_max'].value
# Read arrays.
@ -463,27 +347,15 @@ class WindowedMultipole(EqualityMixin):
out.data = group['data'].value
out.l_value = group['l_value'].value
if out.l_value.shape[0] != out.data.shape[0]:
raise ValueError(err.format('l_value', 'data'))
out.pseudo_k0RS = group['pseudo_K0RS'].value
if out.pseudo_k0RS.shape[0] != out.num_l:
raise ValueError(err.format('pseudo_k0RS', 'l_value'))
out.w_start = group['w_start'].value
out.w_end = group['w_end'].value
if out.w_end.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('w_end', 'w_start'))
out.windows = group['windows'].value
out.broaden_poly = group['broaden_poly'].value.astype(np.bool)
if out.broaden_poly.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('broaden_poly', 'w_start'))
if out.broaden_poly.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('broaden_poly', 'windows'))
out.curvefit = group['curvefit'].value
if out.curvefit.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('curvefit', 'w_start'))
if out.curvefit.shape[0] != out.windows.shape[0]:
raise ValueError(err.format('curvefit', 'windows'))
# _broaden_wmp_polynomials assumes the curve fit has at least 3 terms.
if out.fit_order < 2:
@ -493,7 +365,7 @@ class WindowedMultipole(EqualityMixin):
return out
def _evaluate(self, E, T):
"""Compute total, absorption, and fission cross sections.
"""Compute scattering, absorption, and fission cross sections.
Parameters
----------
@ -510,8 +382,8 @@ class WindowedMultipole(EqualityMixin):
"""
if E < self.start_E: return (0, 0, 0)
if E > self.end_E: return (0, 0, 0)
if E < self.E_min: return (0, 0, 0)
if E > self.E_max: return (0, 0, 0)
# ======================================================================
# Bookkeeping
@ -525,34 +397,12 @@ class WindowedMultipole(EqualityMixin):
# the 1-based vs. 0-based indexing. Similarly startw needs to be
# decreased by 1. endw does not need to be decreased because
# range(startw, endw) does not include endw.
i_window = int(np.floor((sqrtE - sqrt(self.start_E)) / self.spacing))
startw = self.w_start[i_window] - 1
endw = self.w_end[i_window]
# Fill in factors. Because of the unique interference dips in scatering
# resonances, the total cross section has a special "factor" that does
# not appear in the absorption and fission equations.
if startw <= endw:
twophi = np.zeros(self.num_l, dtype=np.float)
sig_t_factor = np.zeros(self.num_l, dtype=np.cfloat)
for iL in range(self.num_l):
twophi[iL] = self.pseudo_k0RS[iL] * sqrtE
if iL == 1:
twophi[iL] = twophi[iL] - np.arctan(twophi[iL])
elif iL == 2:
arg = 3.0 * twophi[iL] / (3.0 - twophi[iL]**2)
twophi[iL] = twophi[iL] - np.arctan(arg)
elif iL == 3:
arg = (twophi[iL] * (15.0 - twophi[iL]**2)
/ (15.0 - 6.0 * twophi[iL]**2))
twophi[iL] = twophi[iL] - np.arctan(arg)
twophi = 2.0 * twophi
sig_t_factor = np.cos(twophi) - 1j*np.sin(twophi)
i_window = int(np.floor((sqrtE - sqrt(self.E_min)) / self.spacing))
startw = self.windows[i_window, 0] - 1
endw = self.windows[i_window, 1]
# Initialize the ouptut cross sections.
sig_t = 0.0
sig_s = 0.0
sig_a = 0.0
sig_f = 0.0
@ -565,7 +415,7 @@ class WindowedMultipole(EqualityMixin):
broadened_polynomials = _broaden_wmp_polynomials(E, dopp,
self.fit_order + 1)
for i_poly in range(self.fit_order+1):
sig_t += (self.curvefit[i_window, i_poly, _FIT_T]
sig_s += (self.curvefit[i_window, i_poly, _FIT_S]
* broadened_polynomials[i_poly])
sig_a += (self.curvefit[i_window, i_poly, _FIT_A]
* broadened_polynomials[i_poly])
@ -575,7 +425,7 @@ class WindowedMultipole(EqualityMixin):
else:
temp = invE
for i_poly in range(self.fit_order+1):
sig_t += self.curvefit[i_window, i_poly, _FIT_T] * temp
sig_s += self.curvefit[i_window, i_poly, _FIT_S] * temp
sig_a += self.curvefit[i_window, i_poly, _FIT_A] * temp
if self.fissionable:
sig_f += self.curvefit[i_window, i_poly, _FIT_F] * temp
@ -589,22 +439,10 @@ class WindowedMultipole(EqualityMixin):
for i_pole in range(startw, endw):
psi_chi = -1j / (self.data[i_pole, _MP_EA] - sqrtE)
c_temp = psi_chi / E
if self.formalism == 'MLBW':
sig_t += ((self.data[i_pole, _MLBW_RT] * c_temp *
sig_t_factor[self.l_value[i_pole]-1]).real
+ (self.data[i_pole, _MLBW_RX] * c_temp).real)
sig_a += (self.data[i_pole, _MLBW_RA] * c_temp).real
if self.fissionable:
sig_f += (self.data[i_pole, _MLBW_RF] * c_temp).real
elif self.formalism == 'RM':
sig_t += (self.data[i_pole, _RM_RT] * c_temp *
sig_t_factor[self.l_value[i_pole]-1]).real
sig_a += (self.data[i_pole, _RM_RA] * c_temp).real
if self.fissionable:
sig_f += (self.data[i_pole, _RM_RF] * c_temp).real
else:
raise ValueError('Unrecognized/Unsupported R-matrix'
' formalism')
sig_s += (self.data[i_pole, _MP_RS] * c_temp).real
sig_a += (self.data[i_pole, _MP_RA] * c_temp).real
if self.fissionable:
sig_f += (self.data[i_pole, _MP_RF] * c_temp).real
else:
# At temperature, use Faddeeva function-based form.
@ -612,27 +450,15 @@ class WindowedMultipole(EqualityMixin):
for i_pole in range(startw, endw):
Z = (sqrtE - self.data[i_pole, _MP_EA]) * dopp
w_val = _faddeeva(Z) * dopp * invE * sqrt(pi)
if self.formalism == 'MLBW':
sig_t += ((self.data[i_pole, _MLBW_RT] *
sig_t_factor[self.l_value[i_pole]-1] +
self.data[i_pole, _MLBW_RX]) * w_val).real
sig_a += (self.data[i_pole, _MLBW_RA] * w_val).real
if self.fissionable:
sig_f += (self.data[i_pole, _MLBW_RF] * w_val).real
elif self.formalism == 'RM':
sig_t += (self.data[i_pole, _RM_RT] * w_val *
sig_t_factor[self.l_value[i_pole]-1]).real
sig_a += (self.data[i_pole, _RM_RA] * w_val).real
if self.fissionable:
sig_f += (self.data[i_pole, _RM_RF] * w_val).real
else:
raise ValueError('Unrecognized/Unsupported R-matrix'
' formalism')
sig_s += (self.data[i_pole, _MP_RS] * w_val).real
sig_a += (self.data[i_pole, _MP_RA] * w_val).real
if self.fissionable:
sig_f += (self.data[i_pole, _MP_RF] * w_val).real
return sig_t, sig_a, sig_f
return sig_s, sig_a, sig_f
def __call__(self, E, T):
"""Compute total, absorption, and fission cross sections.
"""Compute scattering, absorption, and fission cross sections.
Parameters
----------
@ -674,24 +500,14 @@ class WindowedMultipole(EqualityMixin):
g = f.create_group('nuclide')
# Write scalars.
if self.formalism == 'MLBW':
g.create_dataset('formalism',
data=np.array(_FORM_MLBW, dtype=np.int32))
else:
# Assume RM.
g.create_dataset('formalism',
data=np.array(_FORM_RM, dtype=np.int32))
g.create_dataset('spacing', data=np.array(self.spacing))
g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR))
g.create_dataset('start_E', data=np.array(self.start_E))
g.create_dataset('end_E', data=np.array(self.end_E))
g.create_dataset('E_min', data=np.array(self.E_min))
g.create_dataset('E_max', data=np.array(self.E_max))
# Write arrays.
g.create_dataset('data', data=self.data)
g.create_dataset('l_value', data=self.l_value)
g.create_dataset('pseudo_K0RS', data=self.pseudo_k0RS)
g.create_dataset('w_start', data=self.w_start)
g.create_dataset('w_end', data=self.w_end)
g.create_dataset('windows', data=self.windows)
g.create_dataset('broaden_poly',
data=self.broaden_poly.astype(np.int8))
g.create_dataset('curvefit', data=self.curvefit)

View file

@ -23,7 +23,7 @@ except ImportError:
import scipy.sparse as sp
import openmc.data
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from .nuclide import Nuclide, DecayTuple, ReactionTuple
@ -356,7 +356,7 @@ class Chain(object):
if _have_lxml:
tree.write(str(filename), encoding='utf-8', pretty_print=True)
else:
clean_xml_indentation(root_elem)
clean_indentation(root_elem)
tree.write(str(filename), encoding='utf-8')
def form_matrix(self, rates):

View file

@ -1,10 +1,13 @@
from collections import OrderedDict
from collections import OrderedDict, defaultdict
from collections.abc import Iterable
from copy import deepcopy
from pathlib import Path
from xml.etree import ElementTree as ET
import numpy as np
import openmc
from openmc.clean_xml import clean_xml_indentation
import openmc._xml as xml
from openmc.checkvalue import check_type
@ -92,12 +95,104 @@ class Geometry(object):
x.tag, int(x.get('id'))))
# Clean the indentation in the file to be user-readable
clean_xml_indentation(root_element)
xml.clean_indentation(root_element)
# Write the XML Tree to the geometry.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml(cls, path='geometry.xml', materials=None):
"""Generate geometry from XML file
Parameters
----------
path : str, optional
Path to geometry XML file
materials : openmc.Materials or None
Materials used to assign to cells. If None, an attempt is made to
generate it from the materials.xml file.
Returns
-------
openmc.Geometry
Geometry object
"""
# Helper function for keeping a cache of Universe instances
universes = {}
def get_universe(univ_id):
if univ_id not in universes:
univ = openmc.Universe(univ_id)
universes[univ_id] = univ
return universes[univ_id]
tree = ET.parse(path)
root = tree.getroot()
# Get surfaces
surfaces = {}
periodic = {}
for surface in root.findall('surface'):
s = openmc.Surface.from_xml_element(surface)
surfaces[s.id] = s
# Check for periodic surface
other_id = xml.get_text(surface, 'periodic_surface_id')
if other_id is not None:
periodic[s.id] = int(other_id)
# Apply periodic surfaces
for s1, s2 in periodic.items():
surfaces[s1].periodic_surface = surfaces[s2]
# Dictionary that maps each universe to a list of cells/lattices that
# contain it (needed to determine which universe is the root)
child_of = defaultdict(list)
for elem in root.findall('lattice'):
lat = openmc.RectLattice.from_xml_element(elem, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
for u in lat.universes.ravel():
child_of[u].append(lat)
for elem in root.findall('hex_lattice'):
lat = openmc.HexLattice.from_xml_element(elem, get_universe)
universes[lat.id] = lat
if lat.outer is not None:
child_of[lat.outer].append(lat)
if lat.ndim == 2:
for ring in lat.universes:
for u in ring:
child_of[u].append(lat)
else:
for axial_slice in lat.universes:
for ring in axial_slice:
for u in ring:
child_of[u].append(lat)
# Create dictionary to easily look up materials
if materials is None:
filename = Path(path).parent / 'materials.xml'
materials = openmc.Materials.from_xml(str(filename))
mats = {str(m.id): m for m in materials}
mats['void'] = None
for elem in root.findall('cell'):
c = openmc.Cell.from_xml_element(elem, surfaces, mats, get_universe)
if c.fill_type in ('universe', 'lattice'):
child_of[c.fill].append(c)
# Determine which universe is the root by finding one which is not a
# child of any other object
for u in universes.values():
if not child_of[u]:
return cls(u)
else:
raise ValueError('Error determining root universe.')
def find(self, point):
"""Find cells/universes/lattices which contain a given point

View file

@ -10,6 +10,7 @@ import numpy as np
import openmc.checkvalue as cv
import openmc
from openmc._xml import get_text
from openmc.mixin import IDManagerMixin
@ -768,6 +769,42 @@ class RectLattice(Lattice):
# Append the XML subelement for this Lattice to the XML element
xml_element.append(lattice_subelement)
@classmethod
def from_xml_element(cls, elem, get_universe):
"""Generate rectangular lattice from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
`<lattice>` element
get_universe : function
Function returning universe (defined in
:meth:`openmc.Geometry.from_xml`)
Returns
-------
RectLattice
Rectangular lattice
"""
lat_id = int(get_text(elem, 'id'))
name = get_text(elem, 'name')
lat = cls(lat_id, name)
lat.lower_left = [float(i) for i in get_text(elem, 'lower_left').split()]
lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()]
outer = get_text(elem, 'outer')
if outer is not None:
lat.outer = get_universe(int(outer))
# Get array of universes
dimension = get_text(elem, 'dimension').split()
shape = np.array(dimension, dtype=int)[::-1]
uarray = np.array([get_universe(int(i)) for i in
get_text(elem, 'universes').split()])
uarray.shape = shape
lat.universes = uarray
return lat
class HexLattice(Lattice):
r"""A lattice consisting of hexagonal prisms.
@ -1207,6 +1244,78 @@ class HexLattice(Lattice):
# Append the XML subelement for this Lattice to the XML element
xml_element.append(lattice_subelement)
@classmethod
def from_xml_element(cls, elem, get_universe):
"""Generate hexagonal lattice from XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
`<hex_lattice>` element
get_universe : function
Function returning universe (defined in
:meth:`openmc.Geometry.from_xml`)
Returns
-------
HexLattice
Hexagonal lattice
"""
lat_id = int(get_text(elem, 'id'))
name = get_text(elem, 'name')
lat = cls(lat_id, name)
lat.center = [float(i) for i in get_text(elem, 'center').split()]
lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()]
outer = get_text(elem, 'outer')
if outer is not None:
lat.outer = get_universe(int(outer))
# Get nested lists of universes
lat._num_rings = n_rings = int(get_text(elem, 'n_rings'))
lat._num_axial = n_axial = int(get_text(elem, 'n_axial', 1))
# Create empty nested lists for one axial level
univs = [[None for _ in range(max(6*(n_rings - 1 - r), 1))]
for r in range(n_rings)]
if n_axial > 1:
univs = [deepcopy(univs) for i in range(n_axial)]
# Get flat array of universes numbers
uarray = np.array([get_universe(int(i)) for i in
get_text(elem, 'universes').split()])
# Fill nested lists
j = 0
for z in range(n_axial):
# Get list for a single axial level
axial_level = univs[z] if n_axial > 1 else univs
# Start iterating from top
x, alpha = 0, n_rings - 1
while True:
# Set entry in list based on (x,alpha,z) coordinates
_, i_ring, i_within = lat.get_universe_index((x, alpha, z))
axial_level[i_ring][i_within] = uarray[j]
# Move to the right
x += 2
alpha -= 1
if not lat.is_valid_index((x, alpha, z)):
# Move down in y direction
alpha += x - 1
x = 1 - x
if not lat.is_valid_index((x, alpha, z)):
# Move to the right
x += 2
alpha -= 1
if not lat.is_valid_index((x, alpha, z)):
# Reached the bottom
break
j += 1
lat.universes = univs
return lat
def _repr_axial_slice(self, universes):
"""Return string representation for the given 2D group of universes.

View file

@ -9,7 +9,7 @@ import numpy as np
import openmc
import openmc.data
import openmc.checkvalue as cv
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from .mixin import IDManagerMixin
@ -284,6 +284,8 @@ class Material(IDManagerMixin):
# Create the Material
material = cls(mat_id, name)
material.depletable = bool(group.attrs['depletable'])
if 'volume' in group.attrs:
material.volume = group.attrs['volume']
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in group:
@ -833,6 +835,9 @@ class Material(IDManagerMixin):
if self._depletable:
element.set("depletable", "true")
if self._volume:
element.set("volume", str(self._volume))
# Create temperature XML subelement
if self.temperature is not None:
subelement = ET.SubElement(element, "temperature")
@ -907,6 +912,56 @@ class Material(IDManagerMixin):
return element
@classmethod
def from_xml_element(cls, elem):
"""Generate material from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.Material
Material generated from XML element
"""
mat_id = int(elem.get('id'))
mat = cls(mat_id)
mat.name = elem.get('name')
mat.temperature = elem.get('temperature')
mat.depletable = bool(elem.get('depletable'))
# Get each nuclide
for nuclide in elem.findall('nuclide'):
name = nuclide.attrib['name']
if 'ao' in nuclide.attrib:
mat.add_nuclide(name, float(nuclide.attrib['ao']))
elif 'wo' in nuclide.attrib:
mat.add_nuclide(name, float(nuclide.attrib['wo']), 'wo')
# Get each S(a,b) table
for sab in elem.findall('sab'):
fraction = float(sab.get('fraction', 1.0))
mat.add_s_alpha_beta(sab.get('name'), fraction)
# Get total material density
density = elem.find('density')
units = density.get('units')
if units == 'sum':
mat.set_density(units)
else:
value = float(density.get('value'))
mat.set_density(units, value)
# Check for isotropic scattering nuclides
isotropic = elem.find('isotropic')
if isotropic is not None:
mat.isotropic = isotropic.text.split()
return mat
class Materials(cv.CheckedList):
"""Collection of Materials used for an OpenMC simulation.
@ -1026,8 +1081,41 @@ class Materials(cv.CheckedList):
self._create_material_subelements(root_element)
# Clean the indentation in the file to be user-readable
clean_xml_indentation(root_element)
clean_indentation(root_element)
# Write the XML Tree to the materials.xml file
tree = ET.ElementTree(root_element)
tree.write(path, xml_declaration=True, encoding='utf-8')
@classmethod
def from_xml(cls, path='materials.xml'):
"""Generate materials collection from XML file
Parameters
----------
path : str, optional
Path to materials XML file
Returns
-------
openmc.Materials
Materials collection
"""
tree = ET.parse(path)
root = tree.getroot()
# Generate each material
materials = cls()
for material in root.findall('material'):
materials.append(Material.from_xml_element(material))
# Check for cross sections settings
xs = tree.find('cross_sections')
if xs is not None:
materials.cross_sections = xs.text
mpl = tree.find('multipole_library')
if mpl is not None:
materials.multipole_library = mpl.text
return materials

View file

@ -9,7 +9,7 @@ import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from openmc.mixin import IDManagerMixin
@ -816,7 +816,7 @@ class Plots(cv.CheckedList):
self._create_plot_subelements()
# Clean the indentation in the file to be user-readable
clean_xml_indentation(self._plots_file)
clean_indentation(self._plots_file)
# Write the XML Tree to the plots.xml file
tree = ET.ElementTree(self._plots_file)

77
openmc/polynomial.py Normal file
View file

@ -0,0 +1,77 @@
import numpy as np
import openmc
import openmc.capi as capi
def legendre_from_expcoef(coef, domain= (-1,1)):
"""Return a Legendre series object based on expansion coefficients.
Given a list of coefficients from FET tally and a array of down, return
the numpy Legendre object.
Parameters
----------
coef : Iterable of float
A list of coefficients of each term in Legendre polynomials
domain : (2,) List of float
Domain of the Legendre polynomial
Returns
-------
numpy.polynomial.Legendre
A numpy Legendre series class
"""
n = np.arange(len(coef))
c = (2*n + 1) * np.asarray(coef) / (domain[1] - domain[0])
return np.polynomial.Legendre(c, domain)
class Polynomial(object):
"""Abstract Polynomial Class for creating polynomials.
"""
def __init__(self, coef):
self.coef = np.asarray(coef)
class ZernikeRadial(Polynomial):
"""Create radial only Zernike polynomials given coefficients and domain.
The radial only Zernike polynomials are defined as in
:class:`ZernikeRadialFilter`.
Parameters
----------
coef : Iterable of float
A list of coefficients of each term in radial only Zernike polynomials
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
r : float
Position to be evaluated, normalized on radius [0,1]
Attributes
----------
order : int
The maximum (even) order of Zernike polynomials.
radius : float
Domain of Zernike polynomials to be applied on. Default is 1.
norm_coef : iterable of float
The list of coefficients of each term in the polynomials after
normailization.
"""
def __init__(self, coef, radius=1):
super().__init__(coef)
self._order = 2 * (len(self.coef) - 1)
self.radius = radius
norm_vec = (2 * np.arange(len(self.coef)) + 1) / (np.pi * radius**2)
self._norm_coef = norm_vec * self.coef
@property
def order(self):
return self._order
def __call__(self, r):
zn_rad = capi.calc_zn_rad(self.order, r)
return np.sum(self._norm_coef * zn_rad)

View file

@ -6,7 +6,7 @@ import sys
import numpy as np
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
import openmc.checkvalue as cv
from openmc import VolumeCalculation, Source, Mesh
@ -994,7 +994,7 @@ class Settings(object):
self._create_log_grid_bins_subelement(root_element)
# Clean the indentation in the file to be user-readable
clean_xml_indentation(root_element)
clean_indentation(root_element)
# Write the XML Tree to the settings.xml file
tree = ET.ElementTree(root_element)

View file

@ -59,7 +59,6 @@ class Surface(IDManagerMixin):
def __init__(self, surface_id=None, boundary_type='transmission', name=''):
self.id = surface_id
self.name = name
self._type = ''
self.boundary_type = boundary_type
# A dictionary of the quadratic surface coefficients
@ -67,10 +66,6 @@ class Surface(IDManagerMixin):
# Value - coefficient value
self._coefficients = {}
# An ordered list of the coefficient names to export to XML in the
# proper order
self._coeff_keys = []
def __neg__(self):
return Halfspace(self, '-')
@ -203,6 +198,49 @@ class Surface(IDManagerMixin):
return element
@staticmethod
def from_xml_element(elem):
"""Generate surface from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.Surface
Instance of a surface subclass
"""
# Determine appropriate class
surf_type = elem.get('type')
surface_classes = {
'plane': Plane,
'x-plane': XPlane,
'y-plane': YPlane,
'z-plane': ZPlane,
'x-cylinder': XCylinder,
'y-cylinder': YCylinder,
'z-cylinder': ZCylinder,
'sphere': Sphere,
'x-cone': XCone,
'y-cone': YCone,
'z-cone': ZCone,
'quadric': Quadric,
}
cls = surface_classes[surf_type]
# Determine ID, boundary type, coefficients
kwargs = {}
kwargs['surface_id'] = int(elem.get('id'))
kwargs['boundary_type'] = elem.get('boundary', 'transmission')
coeffs = [float(x) for x in elem.get('coeffs').split()]
kwargs.update(dict(zip(cls._coeff_keys, coeffs)))
return cls(**kwargs)
@staticmethod
def from_hdf5(group):
"""Create surface from HDF5 group
@ -324,12 +362,12 @@ class Plane(Surface):
"""
_type = 'plane'
_coeff_keys = ('A', 'B', 'C', 'D')
def __init__(self, surface_id=None, boundary_type='transmission',
A=1., B=0., C=0., D=0., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._type = 'plane'
self._coeff_keys = ['A', 'B', 'C', 'D']
self._periodic_surface = None
self.a = A
self.b = B
@ -458,12 +496,12 @@ class XPlane(Plane):
"""
_type = 'x-plane'
_coeff_keys = ('x0',)
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._type = 'x-plane'
self._coeff_keys = ['x0']
self.x0 = x0
@property
@ -563,13 +601,13 @@ class YPlane(Plane):
"""
_type = 'y-plane'
_coeff_keys = ('y0',)
def __init__(self, surface_id=None, boundary_type='transmission',
y0=0., name=''):
# Initialize YPlane class attributes
super().__init__(surface_id, boundary_type, name=name)
self._type = 'y-plane'
self._coeff_keys = ['y0']
self.y0 = y0
@property
@ -669,13 +707,13 @@ class ZPlane(Plane):
"""
_type = 'z-plane'
_coeff_keys = ('z0',)
def __init__(self, surface_id=None, boundary_type='transmission',
z0=0., name=''):
# Initialize ZPlane class attributes
super().__init__(surface_id, boundary_type, name=name)
self._type = 'z-plane'
self._coeff_keys = ['z0']
self.z0 = z0
@property
@ -774,8 +812,6 @@ class Cylinder(Surface, metaclass=ABCMeta):
def __init__(self, surface_id=None, boundary_type='transmission',
R=1., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['R']
self.r = R
@property
@ -831,12 +867,12 @@ class XCylinder(Cylinder):
"""
_type = 'x-cylinder'
_coeff_keys = ('y0', 'z0', 'R')
def __init__(self, surface_id=None, boundary_type='transmission',
y0=0., z0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'x-cylinder'
self._coeff_keys = ['y0', 'z0', 'R']
self.y0 = y0
self.z0 = z0
@ -953,12 +989,12 @@ class YCylinder(Cylinder):
"""
_type = 'y-cylinder'
_coeff_keys = ('x0', 'z0', 'R')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., z0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'y-cylinder'
self._coeff_keys = ['x0', 'z0', 'R']
self.x0 = x0
self.z0 = z0
@ -1075,12 +1111,12 @@ class ZCylinder(Cylinder):
"""
_type = 'z-cylinder'
_coeff_keys = ('x0', 'y0', 'R')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, R, name=name)
self._type = 'z-cylinder'
self._coeff_keys = ['x0', 'y0', 'R']
self.x0 = x0
self.y0 = y0
@ -1201,12 +1237,12 @@ class Sphere(Surface):
"""
_type = 'sphere'
_coeff_keys = ('x0', 'y0', 'z0', 'R')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R=1., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._type = 'sphere'
self._coeff_keys = ['x0', 'y0', 'z0', 'R']
self.x0 = x0
self.y0 = y0
self.z0 = z0
@ -1348,11 +1384,12 @@ class Cone(Surface, metaclass=ABCMeta):
Type of the surface
"""
_coeff_keys = ('x0', 'y0', 'z0', 'R2')
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
self.x0 = x0
self.y0 = y0
self.z0 = z0
@ -1443,12 +1480,7 @@ class XCone(Cone):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super().__init__(surface_id, boundary_type, x0, y0,
z0, R2, name=name)
self._type = 'x-cone'
_type = 'x-cone'
def evaluate(self, point):
"""Evaluate the surface equation at a given point.
@ -1519,12 +1551,7 @@ class YCone(Cone):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super().__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'y-cone'
_type = 'y-cone'
def evaluate(self, point):
"""Evaluate the surface equation at a given point.
@ -1595,12 +1622,7 @@ class ZCone(Cone):
"""
def __init__(self, surface_id=None, boundary_type='transmission',
x0=0., y0=0., z0=0., R2=1., name=''):
super().__init__(surface_id, boundary_type, x0, y0, z0,
R2, name=name)
self._type = 'z-cone'
_type = 'z-cone'
def evaluate(self, point):
"""Evaluate the surface equation at a given point.
@ -1659,13 +1681,13 @@ class Quadric(Surface):
"""
_type = 'quadric'
_coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k')
def __init__(self, surface_id=None, boundary_type='transmission',
a=0., b=0., c=0., d=0., e=0., f=0., g=0.,
h=0., j=0., k=0., name=''):
super().__init__(surface_id, boundary_type, name=name)
self._type = 'quadric'
self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']
self.a = a
self.b = b
self.c = c

View file

@ -15,7 +15,7 @@ import h5py
import openmc
import openmc.checkvalue as cv
from openmc.clean_xml import clean_xml_indentation
from openmc._xml import clean_indentation
from .mixin import IDManagerMixin
@ -3187,7 +3187,7 @@ class Tallies(cv.CheckedList):
self._create_derivative_subelements(root_element)
# Clean the indentation in the file to be user-readable
clean_xml_indentation(root_element)
clean_indentation(root_element)
# Write the XML Tree to the tallies.xml file
tree = ET.ElementTree(root_element)

2
pyproject.toml Normal file
View file

@ -0,0 +1,2 @@
[build-system]
requires = ["setuptools", "wheel", "numpy", "cython"]

View file

@ -29,9 +29,9 @@ parser.add_argument('-b', '--batch', action='store_true',
args = parser.parse_args()
baseUrl = 'https://github.com/smharper/windowed_multipole_library/blob/master/'
files = ['multipole_lib.tar.gz?raw=true']
checksums = ['3985aea96f7162a9419c7ed8352e6abb']
baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/'
files = ['WMP_Library_v1.0.tar.gz']
checksums = ['22cb675734cfccb278dffd40dcfbf26a']
block_size = 16384
# ==============================================================================
@ -101,12 +101,7 @@ for f in files:
# Extract files
with tarfile.open(fname, 'r') as tgz:
print('Extracting {0}...'.format(fname))
tgz.extractall(path='wmp/')
# Move data files down one level
for filename in glob.glob('wmp/multipole_lib/*'):
shutil.move(filename, 'wmp/')
os.rmdir('wmp/multipole_lib')
tgz.extractall(path='')
# ==============================================================================
# PROMPT USER TO DELETE .TAR.GZ FILES

View file

@ -1,38 +0,0 @@
module angleenergy_header
use hdf5_interface, only: HID_T
!===============================================================================
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
! distribution that is a function of incoming energy. Each derived type must
! implement a sample() subroutine that returns an outgoing energy and scattering
! cosine given an incoming energy.
!===============================================================================
type, abstract :: AngleEnergy
contains
procedure(angleenergy_sample_), deferred :: sample
procedure(angleenergy_from_hdf5_), deferred :: from_hdf5
end type AngleEnergy
abstract interface
subroutine angleenergy_sample_(this, E_in, E_out, mu)
import AngleEnergy
class(AngleEnergy), intent(in) :: this
real(8), intent(in) :: E_in
real(8), intent(out) :: E_out
real(8), intent(out) :: mu
end subroutine angleenergy_sample_
subroutine angleenergy_from_hdf5_(this, group_id)
import AngleEnergy, HID_T
class(AngleEnergy), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
end subroutine angleenergy_from_hdf5_
end interface
type :: AngleEnergyContainer
class(AngleEnergy), allocatable :: obj
end type AngleEnergyContainer
end module angleenergy_header

View file

@ -25,7 +25,7 @@ module openmc_api
use tally_header
use tally_filter_header
use tally_filter
use tally, only: openmc_tally_set_type
use tally, only: openmc_tally_allocate
use simulation
use string, only: to_f_string
use timer_header
@ -50,7 +50,7 @@ module openmc_api
public :: openmc_filter_set_id
public :: openmc_filter_set_type
public :: openmc_finalize
public :: openmc_find
public :: openmc_find_cell
public :: openmc_get_cell_index
public :: openmc_get_keff
public :: openmc_get_filter_index
@ -59,6 +59,7 @@ module openmc_api
public :: openmc_get_nuclide_index
public :: openmc_get_seed
public :: openmc_get_tally_index
public :: openmc_get_tally_next_id
public :: openmc_global_tallies
public :: openmc_hard_reset
public :: openmc_init_f
@ -82,12 +83,16 @@ module openmc_api
public :: openmc_simulation_init
public :: openmc_source_bank
public :: openmc_source_set_strength
public :: openmc_tally_allocate
public :: openmc_tally_get_estimator
public :: openmc_tally_get_id
public :: openmc_tally_get_filters
public :: openmc_tally_get_n_realizations
public :: openmc_tally_get_nuclides
public :: openmc_tally_get_scores
public :: openmc_tally_get_type
public :: openmc_tally_results
public :: openmc_tally_set_estimator
public :: openmc_tally_set_filters
public :: openmc_tally_set_id
public :: openmc_tally_set_nuclides
@ -182,13 +187,12 @@ contains
end function openmc_finalize
!===============================================================================
! OPENMC_FIND determines the ID or a cell or material at a given point in space
! OPENMC_FIND_CELL determines what cell contains a given point in space
!===============================================================================
function openmc_find(xyz, rtype, id, instance) result(err) bind(C)
function openmc_find_cell(xyz, index, instance) result(err) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3) ! Cartesian point
integer(C_INT), intent(in), value :: rtype ! 1 for cell, 2 for material
integer(C_INT32_T), intent(out) :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT32_T), intent(out) :: instance
integer(C_INT) :: err
@ -200,30 +204,22 @@ contains
p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
call find_cell(p, found)
id = -1
index = -1
instance = -1
err = E_UNASSIGNED
if (found) then
if (rtype == 1) then
id = cells(p % coord(p % n_coord) % cell + 1) % id()
elseif (rtype == 2) then
if (p % material == MATERIAL_VOID) then
id = 0
else
id = materials(p % material) % id()
end if
end if
index = p % coord(p % n_coord) % cell + 1
instance = p % cell_instance - 1
err = 0
else
err = E_GEOMETRY
call set_errmsg("Could not find cell/material at position (" // &
call set_errmsg("Could not find cell at position (" // &
trim(to_str(xyz(1))) // "," // trim(to_str(xyz(2))) // "," // &
trim(to_str(xyz(3))) // ").")
end if
end function openmc_find
end function openmc_find_cell
!===============================================================================
! OPENMC_HARD_RESET reset tallies and timers as well as the pseudorandom

View file

@ -1,18 +1,19 @@
#include "cell.h"
#include "openmc/cell.h"
#include <cmath>
#include <sstream>
#include <string>
#include "error.h"
#include "geometry.h"
#include "hdf5_interface.h"
#include "lattice.h"
#include "material.h"
#include "openmc.h"
#include "settings.h"
#include "surface.h"
#include "xml_interface.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -641,7 +642,7 @@ openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
*n = 1;
}
} else {
strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
@ -665,7 +666,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
//TODO: off-by-one
c.material.push_back(i_mat - 1);
} else {
strcpy(openmc_err_msg, "Index in materials array is out of bounds.");
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
@ -676,7 +677,7 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
c.type = FILL_LATTICE;
}
} else {
strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;

View file

@ -243,7 +243,7 @@ contains
use error, only: fatal_error, warning
use mesh_header, only: RegularMesh, openmc_extend_meshes
use string
use tally, only: openmc_tally_set_type
use tally, only: openmc_tally_allocate
use tally_header, only: openmc_extend_tallies
use tally_filter_header
use tally_filter
@ -262,6 +262,7 @@ contains
integer(C_INT) :: err
integer :: i_filt ! index in filters array
integer :: filt_id
integer :: tally_id
integer :: iarray3(3) ! temp integer array
real(8) :: rarray3(3) ! temp double array
real(C_DOUBLE), allocatable :: energies(:)
@ -433,16 +434,13 @@ contains
! Begin loop around tallies
do i = 1, size(cmfd_tallies)
! Allocate tally
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'generic' // C_NULL_CHAR)
err = openmc_tally_allocate(i_start + i - 1, C_CHAR_'generic' // C_NULL_CHAR)
call openmc_get_tally_next_id(tally_id)
err = openmc_tally_set_id(i_start + i - 1, tally_id)
! Point t to tally variable
associate (t => cmfd_tallies(i) % obj)
! Set reset property
if (check_for_node(root, "reset")) then
call get_node_value(root, "reset", t % reset)
end if
! Set the incoming energy mesh filter index in the tally find_filter
! array
n_filter = 1
@ -464,10 +462,10 @@ contains
t % name = "CMFD flux, total"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
t % type = TALLY_VOLUME
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Allocate and set filters
allocate(filter_indices(n_filter))
@ -492,10 +490,10 @@ contains
t % name = "CMFD neutron production"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
t % type = TALLY_VOLUME
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Set the incoming energy mesh filter index in the tally find_filter
! array
@ -527,7 +525,7 @@ contains
t % name = "CMFD surface currents"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Allocate and set filters
allocate(filter_indices(n_filter))
@ -544,17 +542,17 @@ contains
! Set macro bins
t % score_bins(1) = SCORE_CURRENT
t % type = TALLY_MESH_SURFACE
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'mesh-surface' // C_NULL_CHAR)
else if (i == 4) then
! Set name
t % name = "CMFD P1 scatter"
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
t % type = TALLY_VOLUME
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Allocate and set filters
n_filter = 2

View file

@ -25,7 +25,7 @@ module constants
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_VOXEL(2) = [1, 0]
integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0]
character(10), parameter :: VERSION_MULTIPOLE = "v0.2"
character(10), parameter :: VERSION_MULTIPOLE = "v1.0"
! ============================================================================
! ADJUSTABLE PARAMETERS

View file

@ -1,4 +1,4 @@
#include "distribution.h"
#include "openmc/distribution.h"
#include <algorithm> // for copy
#include <cmath> // for sqrt, floor, max
@ -6,10 +6,10 @@
#include <numeric> // for accumulate
#include <string> // for string, stod
#include "error.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "xml_interface.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace openmc {

View file

@ -1,15 +1,16 @@
#include "distribution_angle.h"
#include "openmc/distribution_angle.h"
#include <cmath> // for abs, copysign
#include <vector> // for vector
#include "endf.h"
#include "hdf5_interface.h"
#include "random_lcg.h"
#include "search.h"
#include "xtensor/xarray.hpp"
#include "xtensor/xview.hpp"
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
namespace openmc {
//==============================================================================

View file

@ -1,16 +1,17 @@
#include "distribution_energy.h"
#include "openmc/distribution_energy.h"
#include <algorithm> // for max, min, copy, move
#include <cstddef> // for size_t
#include <iterator> // for back_inserter
#include "endf.h"
#include "hdf5_interface.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "search.h"
#include "xtensor/xview.hpp"
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
namespace openmc {
//==============================================================================

View file

@ -1,11 +1,11 @@
#include "distribution_multi.h"
#include "openmc/distribution_multi.h"
#include <algorithm> // for move
#include <cmath> // for sqrt, sin, cos, max
#include "constants.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "openmc/constants.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
namespace openmc {

View file

@ -1,8 +1,8 @@
#include "distribution_spatial.h"
#include "openmc/distribution_spatial.h"
#include "error.h"
#include "random_lcg.h"
#include "xml_interface.h"
#include "openmc/error.h"
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace openmc {

View file

@ -1,15 +1,16 @@
#include "endf.h"
#include "openmc/endf.h"
#include <algorithm> // for copy
#include <cmath> // for log, exp
#include <iterator> // for back_inserter
#include "constants.h"
#include "hdf5_interface.h"
#include "search.h"
#include "xtensor/xarray.hpp"
#include "xtensor/xview.hpp"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/search.h"
namespace openmc {
//==============================================================================
@ -144,4 +145,35 @@ double Tabulated1D::operator()(double x) const
}
}
//==============================================================================
// CoherentElasticXS implementation
//==============================================================================
CoherentElasticXS::CoherentElasticXS(hid_t dset)
{
// Read 2D array from dataset
xt::xarray<double> arr;
read_dataset(dset, arr);
// Get views for Bragg edges and structure factors
auto E = xt::view(arr, 0);
auto s = xt::view(arr, 1);
// Copy Bragg edges and partial sums of structure factors
std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));
std::copy(s.begin(), s.end(), std::back_inserter(factors_));
}
double CoherentElasticXS::operator()(double E) const
{
if (E < bragg_edges_[0]) {
// If energy is below that of the lowest Bragg peak, the elastic cross
// section will be zero
return 0.0;
} else {
auto i_grid = lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);
return factors_[i_grid] / E;
}
}
} // namespace openmc

View file

@ -1,6 +1,6 @@
#include "finalize.h"
#include "openmc/finalize.h"
#include "message_passing.h"
#include "openmc/message_passing.h"
void openmc_free_bank()
{

View file

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

View file

@ -1,14 +1,14 @@
#include "geometry.h"
#include "openmc/geometry.h"
#include <array>
#include <sstream>
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "lattice.h"
#include "settings.h"
#include "surface.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/lattice.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
namespace openmc {

View file

@ -1,17 +1,17 @@
#include "geometry_aux.h"
#include "openmc/geometry_aux.h"
#include <algorithm> // for std::max
#include <sstream>
#include <unordered_set>
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "geometry.h"
#include "lattice.h"
#include "material.h"
#include "settings.h"
#include "surface.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
namespace openmc {
@ -94,9 +94,9 @@ assign_temperatures()
c->sqrtkT.push_back(0);
} else {
if (global_materials[i_mat]->temperature >= 0) {
if (global_materials[i_mat]->temperature_ >= 0) {
// This material has a default temperature; use that value.
auto T = global_materials[i_mat]->temperature;
auto T = global_materials[i_mat]->temperature_;
c->sqrtkT.push_back(std::sqrt(K_BOLTZMANN * T));
} else {
// Use the global default temperature.

View file

@ -1,4 +1,4 @@
#include "hdf5_interface.h"
#include "openmc/hdf5_interface.h"
#include <array>
#include <cstring>
@ -9,9 +9,9 @@
#include "hdf5_hl.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#include "message_passing.h"
#include "openmc/message_passing.h"
#endif
#include "error.h"
#include "openmc/error.h"
namespace openmc {
@ -327,7 +327,7 @@ get_groups(hid_t group_id, char* name[])
}
std::vector<std::string>
group_names(hid_t group_id)
member_names(hid_t group_id, H5O_type_t type)
{
// Determine number of links in the group
H5G_info_t info;
@ -341,7 +341,7 @@ group_names(hid_t group_id)
// 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;
if (oinfo.type != type) continue;
// Get size of name
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
@ -356,6 +356,17 @@ group_names(hid_t group_id)
return names;
}
std::vector<std::string>
group_names(hid_t group_id)
{
return member_names(group_id, H5O_TYPE_GROUP);
}
std::vector<std::string>
dataset_names(hid_t group_id)
{
return member_names(group_id, H5O_TYPE_DATASET);
}
bool
object_exists(hid_t object_id, const char* name)

View file

@ -1,4 +1,4 @@
#include "initialize.h"
#include "openmc/initialize.h"
#include <cstddef>
#include <cstring>
@ -6,15 +6,16 @@
#include <sstream>
#include <string>
#include "error.h"
#include "hdf5_interface.h"
#include "message_passing.h"
#include "openmc.h"
#include "settings.h"
#ifdef _OPENMP
#include "omp.h"
#endif
#include "openmc/capi.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
// data/functions from Fortran side
extern "C" void print_usage();
extern "C" void print_version();

View file

@ -1398,9 +1398,7 @@ contains
! Check if material is depletable
if (check_for_node(node_mat, "depletable")) then
call get_node_value(node_mat, "depletable", temp_str)
if (to_lower(temp_str) == "true" .or. temp_str == "1") &
mat % depletable = .true.
call get_node_value(node_mat, "depletable", mat % depletable)
end if
! Copy material name
@ -1784,6 +1782,7 @@ contains
integer :: k ! another loop index
integer :: l ! loop over bins
integer :: filter_id ! user-specified identifier for filter
integer :: tally_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
@ -1982,7 +1981,7 @@ contains
READ_TALLIES: do i = 1, n
! Allocate tally
err = openmc_tally_set_type(i_start + i - 1, &
err = openmc_tally_allocate(i_start + i - 1, &
C_CHAR_'generic' // C_NULL_CHAR)
! Get pointer to tally
@ -1991,19 +1990,15 @@ contains
! Get pointer to tally xml node
node_tal = node_tal_list(i)
! Copy tally id
! Copy and set tally id
if (check_for_node(node_tal, "id")) then
call get_node_value(node_tal, "id", t % id)
call get_node_value(node_tal, "id", tally_id)
err = openmc_tally_set_id(i_start + i - 1, tally_id)
if (err /= 0) call fatal_error(to_f_string(openmc_err_msg))
else
call fatal_error("Must specify id for tally in tally XML file.")
end if
! Check to make sure 'id' hasn't been used
if (tally_dict % has(t % id)) then
call fatal_error("Two or more tallies use the same unique ID: " &
// to_str(t % id))
end if
! Copy tally name
if (check_for_node(node_tal, "name")) &
call get_node_value(node_tal, "name", t % name)
@ -2703,9 +2698,6 @@ contains
end select
end if
! Add tally to dictionary
call tally_dict % set(t % id, i)
end associate
end do READ_TALLIES

View file

@ -1,15 +1,15 @@
#include "lattice.h"
#include "openmc/lattice.h"
#include <cmath>
#include <sstream>
#include <vector>
#include "cell.h"
#include "error.h"
#include "geometry_aux.h"
#include "hdf5_interface.h"
#include "string_utils.h"
#include "xml_interface.h"
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/geometry_aux.h"
#include "openmc/hdf5_interface.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
namespace openmc {

View file

@ -1,8 +1,8 @@
#include "error.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "openmc.h"
#include "openmc/capi.h"
#include "openmc/error.h"
int main(int argc, char* argv[]) {

View file

@ -1,10 +1,10 @@
#include "material.h"
#include "openmc/material.h"
#include <string>
#include <sstream>
#include "error.h"
#include "xml_interface.h"
#include "openmc/error.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -20,16 +20,20 @@ std::unordered_map<int32_t, int32_t> material_map;
// Material implementation
//==============================================================================
Material::Material(pugi::xml_node material_node)
Material::Material(pugi::xml_node node)
{
if (check_for_node(material_node, "id")) {
id = std::stoi(get_node_value(material_node, "id"));
if (check_for_node(node, "id")) {
id = std::stoi(get_node_value(node, "id"));
} else {
fatal_error("Must specify id of material in materials XML file.");
}
if (check_for_node(material_node, "temperature")) {
temperature = std::stod(get_node_value(material_node, "temperature"));
if (check_for_node(node, "temperature")) {
temperature_ = std::stod(get_node_value(node, "temperature"));
}
if (check_for_node(node, "volume")) {
volume_ = std::stod(get_node_value(node, "volume"));
}
}
@ -60,6 +64,48 @@ read_materials(pugi::xml_node* node)
}
}
//==============================================================================
// C API
//==============================================================================
extern "C" int
openmc_material_get_volume(int32_t index, double* volume)
{
if (index >= 1 && index <= global_materials.size()) {
Material* m = global_materials[index - 1];
if (m->volume_ >= 0.0) {
*volume = m->volume_;
return 0;
} else {
std::stringstream msg;
msg << "Volume for material with ID=" << m->id << " not set.";
set_errmsg(msg);
return OPENMC_E_UNASSIGNED;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
extern "C" int
openmc_material_set_volume(int32_t index, double volume)
{
if (index >= 1 && index <= global_materials.size()) {
Material* m = global_materials[index - 1];
if (volume >= 0.0) {
m->volume_ = volume;
return 0;
} else {
set_errmsg("Volume must be non-negative");
return OPENMC_E_INVALID_ARGUMENT;
}
} else {
set_errmsg("Index in materials array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================

View file

@ -24,6 +24,7 @@ module material_header
public :: openmc_material_add_nuclide
public :: openmc_material_get_id
public :: openmc_material_get_densities
public :: openmc_material_get_volume
public :: openmc_material_set_density
public :: openmc_material_set_densities
public :: openmc_material_set_id
@ -51,9 +52,16 @@ module material_header
end subroutine material_set_id_c
subroutine extend_materials_c(n) bind(C)
import C_INT32_t
import C_INT32_T
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_materials_c
function openmc_material_get_volume(index, volume) result(err) bind(C)
import C_INT32_T, C_DOUBLE, C_INT
integer(C_INT32_T), value :: index
real(C_DOUBLE), intent(out) :: volume
integer(C_INT) :: err
end function openmc_material_get_volume
end interface
!===============================================================================
@ -215,7 +223,7 @@ contains
found = .false.
associate (sab => sab_tables(this % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(this % nuclide)
if (any(sab % nuclides == nuclides(this % nuclide(j)) % name)) then
if (sab % has_nuclide(nuclides(this % nuclide(j)) % name)) then
call i_sab_tables % push_back(this % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
call sab_fracs % push_back(this % sab_fracs(k))
@ -369,7 +377,7 @@ contains
! If particle energy is greater than the highest energy for the
! S(a,b) table, then don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
if (p % E > sab_tables(i_sab) % threshold()) then
i_sab = 0
end if

View file

@ -1,4 +1,4 @@
#include "math_functions.h"
#include "openmc/math_functions.h"
namespace openmc {
@ -608,7 +608,7 @@ void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
double k2 = 2 * p * (p - 1) * (p - 2);
double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);
double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;
zn_rad[index] =
zn_rad[index] =
((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1;
}
}

View file

@ -1,4 +1,4 @@
#include "message_passing.h"
#include "openmc/message_passing.h"
namespace openmc {
namespace mpi {

View file

@ -1,17 +1,18 @@
#include "openmc/mgxs.h"
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <valarray>
#ifdef _OPENMP
# include <omp.h>
#endif
#ifdef _OPENMP
#include <omp.h>
#endif
#include "error.h"
#include "math_functions.h"
#include "random_lcg.h"
#include "string_functions.h"
#include "mgxs.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/string_functions.h"
namespace openmc {
@ -709,4 +710,4 @@ Mgxs::set_angle_index(const double uvw[3])
}
}
} // namespace openmc
} // namespace openmc

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