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

Merge pull request #1065 from nelsonag/xtensor_mg

Browse source 
Use the xtensor library in the multi-group solver
This commit is contained in:
Paul Romano 2018-09-08 16:09:38 -05:00 • committed by GitHub
commit 74b7e53bbb
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
36 changed files with 1350 additions and 1470 deletions
Download patch file Download diff file Expand all files Collapse all files

13
include/openmc/constants.h
View file

@ -11,16 +11,9 @@
namespace openmc {
// TODO: Replace with xtensor/other library?
typedef std::vector<double> double_1dvec;
typedef std::vector<std::vector<double> > double_2dvec;
typedef std::vector<std::vector<std::vector<double> > > double_3dvec;
typedef std::vector<std::vector<std::vector<std::vector<double> > > > double_4dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > double_5dvec;
typedef std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > > > double_6dvec;
typedef std::vector<int> int_1dvec;
typedef std::vector<std::vector<int> > int_2dvec;
typedef std::vector<std::vector<std::vector<int> > > int_3dvec;
using double_2dvec = std::vector<std::vector<double>>;
using double_3dvec = std::vector<std::vector<std::vector<double>>>;
using double_4dvec = std::vector<std::vector<std::vector<std::vector<double>>>>;
// ============================================================================
// VERSIONING NUMBERS

70
include/openmc/hdf5_interface.h
View file

@ -14,6 +14,7 @@
#include "xtensor/xarray.hpp"
#include "openmc/position.h"
#include "openmc/error.h"
namespace openmc {
@ -46,39 +47,6 @@ hid_t file_open(const std::string& filename, char mode, bool parallel=false);
void write_string(hid_t group_id, const char* name, const std::string& buffer,
bool indep);
void
read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<int> > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<double> > > >& result,
bool must_have = false);
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
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);
@ -236,7 +204,7 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
}
//==============================================================================
// Templates/overloads for read_dataset
// Templates/overloads for read_dataset and related methods
//==============================================================================
template<typename T>
@ -294,6 +262,40 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep
close_dataset(dset);
}
template <typename T, std::size_t N>
void read_dataset_as_shape(hid_t obj_id, const char* name,
xt::xtensor<T, N>& arr, bool indep=false)
{
hid_t dset = open_dataset(obj_id, name);
// Allocate new array to read data into
std::size_t size = 1;
for (const auto x : arr.shape())
size *= x;
T* buffer = new T[size];
// Read data from attribute
read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer, indep);
// Adapt into xarray
arr = xt::adapt(buffer, size, xt::acquire_ownership(), arr.shape());
close_dataset(dset);
}
template <typename T, std::size_t N>
void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
bool must_have=false)
{
if (object_exists(obj_id, name)) {
read_dataset_as_shape(obj_id, name, result, true);
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
//==============================================================================
// Templates/overloads for write_attribute
//==============================================================================

26
include/openmc/mgxs.h
View file

@ -7,6 +7,8 @@
#include <string>
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xsdata.h"
@ -35,7 +37,7 @@ struct CacheData {
class Mgxs {
private:
double_1dvec kTs; // temperature in eV (k * T)
xt::xtensor<double, 1> kTs; // temperature in eV (k * T)
int scatter_format; // flag for if this is legendre, histogram, or tabular
int num_delayed_groups; // number of delayed neutron groups
int num_groups; // number of energy groups
@ -44,8 +46,8 @@ class Mgxs {
bool is_isotropic; // used to skip search for angle indices if isotropic
int n_pol;
int n_azi;
double_1dvec polar;
double_1dvec azimuthal;
std::vector<double> polar;
std::vector<double> azimuthal;
//! \brief Initializes the Mgxs object metadata
//!
@ -62,10 +64,10 @@ class Mgxs {
//! @param in_polar Polar angle grid.
//! @param in_azimuthal Azimuthal angle grid.
void
init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs,
init(const std::string& in_name, double in_awr, const std::vector<double>& in_kTs,
bool in_fissionable, int in_scatter_format, int in_num_groups,
int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal);
const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal);
//! \brief Initializes the Mgxs object metadata from the HDF5 file
//!
@ -80,8 +82,8 @@ class Mgxs {
//! @param method Method of choosing nearest temperatures.
void
metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim,
int in_num_delayed_groups, const std::vector<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim,
int& method);
//! \brief Performs the actual act of combining the microscopic data for a
@ -93,8 +95,8 @@ class Mgxs {
//! corresponds to the temperature of interest.
//! @param this_t The temperature index of the macroscopic object.
void
combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t);
combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scalars,
const std::vector<int>& micro_ts, int this_t);
//! \brief Checks to see if this and that are able to be combined
//!
@ -128,7 +130,7 @@ class Mgxs {
//! provides the number of points to use in the tabular representation.
//! @param method Method of choosing nearest temperatures.
Mgxs(hid_t xs_id, int energy_groups,
int delayed_groups, const double_1dvec& temperature, double tolerance,
int delayed_groups, const std::vector<double>& temperature, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points, int& method);
@ -141,8 +143,8 @@ class Mgxs {
//! @param atom_densities Atom densities of those microscopic quantities.
//! @param tolerance Tolerance of temperature selection method.
//! @param method Method of choosing nearest temperatures.
Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method);
//! \brief Provides a cross section value given certain parameters

77
include/openmc/scattdata.h
View file

@ -6,6 +6,8 @@
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/constants.h"
namespace openmc {
@ -25,23 +27,25 @@ class ScattData {
protected:
//! \brief Initializes the attributes of the base class.
void
base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_energy, const double_2dvec& in_mult);
base_init(int order, const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult);
//! \brief Combines microscopic ScattDatas into a macroscopic one.
void
base_combine(int max_order, const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& sparse_mult, double_3dvec& sparse_scatter);
base_combine(size_t max_order, const std::vector<ScattData*>& those_scatts,
const std::vector<double>& scalars, xt::xtensor<int, 1>& in_gmin,
xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter);
public:
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
int_1dvec gmin; // minimum outgoing group
int_1dvec gmax; // maximum outgoing group
double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}}
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
xt::xtensor<double, 1> gmin; // minimum outgoing group
xt::xtensor<double, 1> gmax; // maximum outgoing group
xt::xtensor<double, 1> scattxs; // Isotropic Sigma_{s,g_{in}}
//! \brief Calculates the value of normalized f(mu).
//!
@ -72,7 +76,7 @@ class ScattData {
//! @param in_mult Input sparse multiplicity matrix
//! @param coeffs Input sparse scattering matrix
virtual void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs) = 0;
//! \brief Combines the microscopic data.
@ -81,7 +85,7 @@ class ScattData {
//! @param scalars Scalars to multiply the microscopic data by.
virtual void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars) = 0;
const std::vector<double>& scalars) = 0;
//! \brief Getter for the dimensionality of the scattering order.
//!
@ -89,7 +93,7 @@ class ScattData {
//! of points, and for Histogram this is the number of bins.
//!
//! @return The order.
virtual int
virtual size_t
get_order() = 0;
//! \brief Builds a dense scattering matrix from the constituent parts
@ -97,8 +101,8 @@ class ScattData {
//! @param max_order If Legendre this is the maximum value of "n" in "Pn"
//! requested; ignored otherwise.
//! @return The dense scattering matrix.
virtual double_3dvec
get_matrix(int max_order) = 0;
virtual xt::xtensor<double, 3>
get_matrix(size_t max_order) = 0;
//! \brief Samples the outgoing energy from the ScattData info.
//!
@ -142,12 +146,12 @@ class ScattDataLegendre: public ScattData {
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const std::vector<double>& scalars);
//! \brief Find the maximal value of the angular distribution to use as a
// bounding box with rejection sampling.
@ -160,11 +164,11 @@ class ScattDataLegendre: public ScattData {
void
sample(int gin, int& gout, double& mu, double& wgt);
int
size_t
get_order() {return dist[0][0].size() - 1;};
double_3dvec
get_matrix(int max_order);
xt::xtensor<double, 3>
get_matrix(size_t max_order);
};
//==============================================================================
@ -176,19 +180,19 @@ class ScattDataHistogram: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the spacing between the mu bin points
double_3dvec fmu; // The angular distribution histogram
xt::xtensor<double, 1> mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the mu spacing
double_3dvec fmu; // The angular distribution histogram
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const std::vector<double>& scalars);
double
calc_f(int gin, int gout, double mu);
@ -196,11 +200,11 @@ class ScattDataHistogram: public ScattData {
void
sample(int gin, int& gout, double& mu, double& wgt);
int
size_t
get_order() {return dist[0][0].size();};
double_3dvec
get_matrix(int max_order);
xt::xtensor<double, 3>
get_matrix(size_t max_order);
};
//==============================================================================
@ -212,9 +216,9 @@ class ScattDataTabular: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu grid points
double dmu; // Quick storage of the spacing between the mu points
double_3dvec fmu; // The angular distribution function
xt::xtensor<double, 1> mu; // Angle distribution mu grid points
double dmu; // Quick storage of the mu spacing
double_3dvec fmu; // The angular distribution function
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
@ -225,12 +229,12 @@ class ScattDataTabular: public ScattData {
public:
void
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
init(const xt::xtensor<int, 1>& in_gmin, const xt::xtensor<int, 1>& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const std::vector<double>& scalars);
double
calc_f(int gin, int gout, double mu);
@ -238,10 +242,11 @@ class ScattDataTabular: public ScattData {
void
sample(int gin, int& gout, double& mu, double& wgt);
int
size_t
get_order() {return dist[0][0].size();};
double_3dvec get_matrix(int max_order);
xt::xtensor<double, 3>
get_matrix(size_t max_order);
};
//==============================================================================

72
include/openmc/xsdata.h
View file

@ -7,6 +7,8 @@
#include <memory>
#include <vector>
#include "xtensor/xtensor.hpp"
#include "openmc/hdf5_interface.h"
#include "openmc/scattdata.h"
@ -22,41 +24,77 @@ class XsData {
private:
//! \brief Reads scattering data from the HDF5 file
void
scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups,
int scatter_format, int final_scatter_format, int order_data,
int max_order, int legendre_to_tabular_points);
//! \brief Reads fission data from the HDF5 file
void
fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
int delayed_groups, bool is_isotropic);
fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups,
size_t delayed_groups, bool is_isotropic);
//! \brief Reads fission data formatted as chi and nu-fission vectors from
// the HDF5 file when beta is provided.
void
fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups, size_t delayed_groups, bool is_isotropic);
//! \brief Reads fission data formatted as chi and nu-fission vectors from
// the HDF5 file when beta is not provided.
void
fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups, size_t delayed_groups);
//! \brief Reads fission data formatted as chi and nu-fission vectors from
// the HDF5 file when no delayed data is provided.
void
fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups);
//! \brief Reads fission data formatted as a nu-fission matrix from
// the HDF5 file when beta is provided.
void
fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups, size_t delayed_groups, bool is_isotropic);
//! \brief Reads fission data formatted as a nu-fission matrix from
// the HDF5 file when beta is not provided.
void
fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups, size_t delayed_groups);
//! \brief Reads fission data formatted as a nu-fission matrix from
// the HDF5 file when no delayed data is provided.
void
fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang,
size_t energy_groups);
public:
// The following quantities have the following dimensions:
// [angle][incoming group]
double_2dvec total;
double_2dvec absorption;
double_2dvec nu_fission;
double_2dvec prompt_nu_fission;
double_2dvec kappa_fission;
double_2dvec fission;
double_2dvec inverse_velocity;
xt::xtensor<double, 2> total;
xt::xtensor<double, 2> absorption;
xt::xtensor<double, 2> nu_fission;
xt::xtensor<double, 2> prompt_nu_fission;
xt::xtensor<double, 2> kappa_fission;
xt::xtensor<double, 2> fission;
xt::xtensor<double, 2> inverse_velocity;
// decay_rate has the following dimensions:
// [angle][delayed group]
double_2dvec decay_rate;
xt::xtensor<double, 2> decay_rate;
// delayed_nu_fission has the following dimensions:
// [angle][incoming group][delayed group]
double_3dvec delayed_nu_fission;
xt::xtensor<double, 3> delayed_nu_fission;
// chi_prompt has the following dimensions:
// [angle][incoming group][outgoing group]
double_3dvec chi_prompt;
xt::xtensor<double, 3> chi_prompt;
// chi_delayed has the following dimensions:
// [angle][incoming group][outgoing group][delayed group]
double_4dvec chi_delayed;
xt::xtensor<double, 4> chi_delayed;
// scatter has the following dimensions: [angle]
std::vector<std::shared_ptr<ScattData> > scatter;
std::vector<std::shared_ptr<ScattData>> scatter;
XsData() = default;
@ -68,7 +106,7 @@ class XsData {
//! @param scatter_format The scattering representation of the file.
//! @param n_pol Number of polar angles.
//! @param n_azi Number of azimuthal angles.
XsData(int num_groups, int num_delayed_groups, bool fissionable,
XsData(size_t num_groups, size_t num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi);
//! \brief Loads the XsData object from the HDF5 file
@ -101,7 +139,7 @@ class XsData {
//! @param micros Microscopic objects to combine.
//! @param scalars Scalars to multiply the microscopic data by.
void
combine(const std::vector<XsData*>& those_xs, const double_1dvec& scalars);
combine(const std::vector<XsData*>& those_xs, const std::vector<double>& scalars);
//! \brief Checks to see if this and that are able to be combined
//!

163
openmc/examples.py
View file

@ -1,3 +1,5 @@
from numbers import Integral
import numpy as np
import openmc
@ -538,20 +540,20 @@ def pwr_assembly():
return model
def slab_mg(reps=None, as_macro=True):
"""Create a one-group, 1D slab model.
def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5'):
"""Create a 1D slab model.
Parameters
----------
reps : list, optional
List of angular representations. Each item corresponds to materials and
dictates the angular representation of the multi-group cross
sections---isotropic ('iso') or angle-dependent ('ang'), and if Legendre
scattering or tabular scattering ('mu') is used. Thus, items can be
'ang', 'ang_mu', 'iso', or 'iso_mu'.
num_regions : int, optional
Number of regions in the problem, each with a unique MGXS dataset.
Defaults to 1.
as_macro : bool, optional
Whether :class:`openmc.Macroscopic` is used
mat_names : Iterable of str, optional
List of the material names to use; defaults to ['mat_1', 'mat_2',...].
mgxslib_name : str, optional
MGXS Library file to use; defaults to '2g.h5'.
Returns
-------
@ -559,71 +561,82 @@ def slab_mg(reps=None, as_macro=True):
One-group, 1D slab model
"""
openmc.check_type('num_regions', num_regions, Integral)
openmc.check_greater_than('num_regions', num_regions, 0)
if mat_names is not None:
openmc.check_length('mat_names', mat_names, num_regions)
openmc.check_iterable_type('mat_names', mat_names, str)
else:
mat_names = []
for i in range(num_regions):
mat_names.append('mat_' + str(i + 1))
# # Make Materials
materials_file = openmc.Materials()
macros = []
mats = []
for i in range(len(mat_names)):
macros.append(openmc.Macroscopic('mat_' + str(i + 1)))
mats.append(openmc.Material(name=mat_names[i]))
mats[-1].set_density('macro', 1.0)
mats[-1].add_macroscopic(macros[-1])
materials_file += mats
materials_file.cross_sections = mgxslib_name
# # Make Geometry
rad_outer = 929.45
# Set a cell boundary to exist for every material above (exclude the 0)
rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:]
# Instantiate Universe
root = openmc.Universe(universe_id=0, name='root universe')
cells = []
surfs = []
surfs.append(openmc.XPlane(x0=0., boundary_type='reflective'))
for r, rad in enumerate(rads):
if r == len(rads) - 1:
surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum'))
else:
surfs.append(openmc.XPlane(x0=rad))
# Instantiate Cells
cells = []
for c in range(len(surfs) - 1):
cells.append(openmc.Cell())
cells[-1].region = (+surfs[c] & -surfs[c + 1])
cells[-1].fill = mats[c]
# Register Cells with Universe
root.add_cells(cells)
# Instantiate a Geometry, register the root Universe, and export to XML
geometry_file = openmc.Geometry(root)
# # Make Settings
# Instantiate a Settings object, set all runtime parameters
settings_file = openmc.Settings()
settings_file.energy_mode = "multi-group"
settings_file.tabular_legendre = {'enable': False}
settings_file.batches = 10
settings_file.inactive = 5
settings_file.particles = 1000
# Build source distribution
INF = 1000.
bounds = [0., -INF, -INF, rads[0], INF, INF]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
settings_file.source = openmc.source.Source(space=uniform_dist)
settings_file.output = {'summary': False}
model = openmc.model.Model()
# Define materials needed for 1D/1G slab problem
mat_names = ['uo2', 'clad', 'lwtr']
mgxs_reps = ['ang', 'ang_mu', 'iso', 'iso_mu']
if reps is None:
reps = mgxs_reps
xs = []
i = 0
for mat in mat_names:
for rep in reps:
i += 1
name = mat + '_' + rep
xs.append(name)
if as_macro:
m = openmc.Material(name=str(i))
m.set_density('macro', 1.)
m.add_macroscopic(name)
else:
m = openmc.Material(name=str(i))
m.set_density('atom/b-cm', 1.)
m.add_nuclide(name, 1.0, 'ao')
model.materials.append(m)
# Define the materials file
model.xs_data = xs
model.materials.cross_sections = "../../1d_mgxs.h5"
# Define surfaces.
# Assembly/Problem Boundary
left = openmc.XPlane(x0=0.0, boundary_type='reflective')
right = openmc.XPlane(x0=10.0, boundary_type='reflective')
bottom = openmc.YPlane(y0=0.0, boundary_type='reflective')
top = openmc.YPlane(y0=10.0, boundary_type='reflective')
# for each material add a plane
planes = [openmc.ZPlane(z0=0.0, boundary_type='reflective')]
dz = round(5. / float(len(model.materials)), 4)
for i in range(len(model.materials) - 1):
planes.append(openmc.ZPlane(z0=dz * float(i + 1)))
planes.append(openmc.ZPlane(z0=5.0, boundary_type='reflective'))
# Define cells for each material
model.geometry.root_universe = openmc.Universe(name='root universe')
xy = +left & -right & +bottom & -top
for i, mat in enumerate(model.materials):
c = openmc.Cell(fill=mat, region=xy & +planes[i] & -planes[i + 1])
model.geometry.root_universe.add_cell(c)
model.settings.batches = 10
model.settings.inactive = 5
model.settings.particles = 100
model.settings.source = openmc.Source(space=openmc.stats.Box(
[0.0, 0.0, 0.0], [10.0, 10.0, 5.]))
model.settings.energy_mode = "multi-group"
plot = openmc.Plot()
plot.filename = 'mat'
plot.origin = (5.0, 5.0, 2.5)
plot.width = (2.5, 2.5)
plot.basis = 'xz'
plot.pixels = (3000, 3000)
plot.color_by = 'material'
model.plots.append(plot)
model.geometry = geometry_file
model.materials = materials_file
model.settings = settings_file
model.xs_data = macros
return model

170
src/hdf5_interface.cpp
View file

@ -5,13 +5,15 @@
#include <sstream>
#include <string>
#include "xtensor/xtensor.hpp"
#include "xtensor/xarray.hpp"
#include "hdf5.h"
#include "hdf5_hl.h"
#ifdef OPENMC_MPI
#include "mpi.h"
#include "openmc/message_passing.h"
#endif
#include "openmc/error.h"
namespace openmc {
@ -532,172 +534,6 @@ read_complex(hid_t obj_id, const char* name, std::complex<double>* buffer, bool
}
void
read_nd_vector(hid_t obj_id, const char* name, std::vector<double>& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
read_double(obj_id, name, result.data(), true);
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<double> >& result, bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
double temp_arr[dim1 * dim2];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
result[i][j] = temp_arr[temp_idx++];
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<int> >& result, bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int temp_arr[dim1 * dim2];
read_int(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
result[i][j] = temp_arr[temp_idx++];
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<double> > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
double temp_arr[dim1 * dim2 * dim3];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<int> > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int temp_arr[dim1 * dim2 * dim3];
read_int(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
result[i][j][k] = temp_arr[temp_idx++];
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<double> > > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
double temp_arr[dim1 * dim2 * dim3 * dim4];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
result[i][j][k][l] = temp_arr[temp_idx++];
}
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_nd_vector(hid_t obj_id, const char* name,
std::vector<std::vector<std::vector<std::vector<std::vector<double> > > > >& result,
bool must_have)
{
if (object_exists(obj_id, name)) {
int dim1 = result.size();
int dim2 = result[0].size();
int dim3 = result[0][0].size();
int dim4 = result[0][0][0].size();
int dim5 = result[0][0][0][0].size();
double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5];
read_double(obj_id, name, temp_arr, true);
int temp_idx = 0;
for (int i = 0; i < dim1; i++) {
for (int j = 0; j < dim2; j++) {
for (int k = 0; k < dim3; k++) {
for (int l = 0; l < dim4; l++) {
for (int m = 0; m < dim5; m++) {
result[i][j][k][l][m] = temp_arr[temp_idx++];
}
}
}
}
}
} else if (must_have) {
fatal_error(std::string("Must provide " + std::string(name) + "!"));
}
}
void
read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results)
{

2
src/input_xml.F90
View file

@ -181,7 +181,7 @@ contains
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures()
! Determine desired txemperatures for each nuclide and S(a,b) table
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
! Check to make sure there are not too many nested coordinate levels in the

172
src/mgxs.cpp
View file

@ -3,12 +3,17 @@
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <valarray>
#include <sstream>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xmath.hpp"
#include "xtensor/xsort.hpp"
#include "xtensor/xadapt.hpp"
#include "xtensor/xview.hpp"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
@ -28,14 +33,15 @@ std::vector<Mgxs> macro_xs;
void
Mgxs::init(const std::string& in_name, double in_awr,
const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format,
const std::vector<double>& in_kTs, bool in_fissionable, int in_scatter_format,
int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal)
const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal)
{
// Set the metadata
name = in_name;
awr = in_awr;
kTs = in_kTs;
//TODO: Remove adapt when in_KTs is an xtensor
kTs = xt::adapt(in_kTs);
fissionable = in_fissionable;
scatter_format = in_scatter_format;
num_groups = in_num_groups;
@ -61,8 +67,8 @@ Mgxs::init(const std::string& in_name, double in_awr,
void
Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method)
int in_num_delayed_groups, const std::vector<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim, int& method)
{
// get name
char char_name[MAX_WORD_LEN];
@ -87,7 +93,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
dset_names[i] = new char[151];
}
get_datasets(kT_group, dset_names);
double_1dvec available_temps(num_temps);
xt::xarray<double> available_temps(num_temps);
for (int i = 0; i < num_temps; i++) {
read_double(kT_group, dset_names[i], &available_temps[i], true);
@ -110,24 +116,20 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
switch(method) {
case TEMPERATURE_NEAREST:
// Find the minimum difference
for (int i = 0; i < temperature.size(); i++) {
std::valarray<double> temp_diff(available_temps.data(),
available_temps.size());
temp_diff = std::abs(temp_diff - temperature[i]);
int i_closest = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// Determine actual temperatures to read
for (const auto& T : temperature) {
auto i_closest = xt::argmin(xt::abs(available_temps - T))[0];
double temp_actual = available_temps[i_closest];
if (std::abs(temp_actual - temperature[i]) < tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(),
std::round(temp_actual)) == temps_to_read.end()) {
if (std::fabs(temp_actual - T) < tolerance) {
if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual))
== temps_to_read.end()) {
temps_to_read.push_back(std::round(temp_actual));
} else {
fatal_error("MGXS Library does not contain cross section for " +
in_name + " at or near " +
std::to_string(std::round(temperature[i])) + " K.");
}
} else {
std::stringstream msg;
msg << "MGXS library does not contain cross sections for "
<< in_name << " at or near " << std::round(T) << " K.";
fatal_error(msg);
}
}
break;
@ -160,7 +162,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
// Get the library's temperatures
int n_temperature = temps_to_read.size();
double_1dvec in_kTs(n_temperature);
std::vector<double> in_kTs(n_temperature);
for (int i = 0; i < n_temperature; i++) {
std::string temp_str(std::to_string(temps_to_read[i]) + "K");
@ -254,12 +256,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
}
// Set the angular bins to use equally-spaced bins
double_1dvec in_polar(in_n_pol);
std::vector<double> in_polar(in_n_pol);
double dangle = PI / in_n_pol;
for (int p = 0; p < in_n_pol; p++) {
in_polar[p] = (p + 0.5) * dangle;
}
double_1dvec in_azimuthal(in_n_azi);
std::vector<double> in_azimuthal(in_n_azi);
dangle = 2. * PI / in_n_azi;
for (int a = 0; a < in_n_azi; a++) {
in_azimuthal[a] = (a + 0.5) * dangle - PI;
@ -274,12 +276,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
//==============================================================================
Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
const double_1dvec& temperature, double tolerance, int max_order,
const std::vector<double>& temperature, double tolerance, int max_order,
bool legendre_to_tabular, int legendre_to_tabular_points, int& method)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
std::vector<int> temps_to_read;
metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
tolerance, temps_to_read, order_data, method);
@ -310,8 +312,8 @@ Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
//==============================================================================
Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method)
{
// Get the minimum data needed to initialize:
@ -328,8 +330,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
int in_num_groups = micros[0]->num_groups;
int in_num_delayed_groups = micros[0]->num_delayed_groups;
bool in_is_isotropic = micros[0]->is_isotropic;
double_1dvec in_polar = micros[0]->polar;
double_1dvec in_azimuthal = micros[0]->azimuthal;
std::vector<double> in_polar = micros[0]->polar;
std::vector<double> in_azimuthal = micros[0]->azimuthal;
init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
@ -345,33 +347,27 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
// Create the list of temperature indices and interpolation factors for
// each microscopic data at the material temperature
int_1dvec micro_t(micros.size(), 0);
double_1dvec micro_t_interp(micros.size(), 0.);
std::vector<int> micro_t(micros.size(), 0);
std::vector<double> micro_t_interp(micros.size(), 0.);
for (int m = 0; m < micros.size(); m++) {
switch(method) {
case TEMPERATURE_NEAREST:
{
// Find the nearest temperature
std::valarray<double> temp_diff(micros[m]->kTs.data(),
micros[m]->kTs.size());
temp_diff = std::abs(temp_diff - temp_desired);
micro_t[m] = std::min_element(std::begin(temp_diff),
std::end(temp_diff)) -
std::begin(temp_diff);
double temp_actual = micros[m]->kTs[micro_t[m]];
micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0];
auto temp_actual = micros[m]->kTs[micro_t[m]];
if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) {
fatal_error("MGXS Library does not contain cross section for " +
name + " at or near " +
std::to_string(std::round(temp_desired / K_BOLTZMANN))
+ " K.");
std::stringstream msg;
msg << "MGXS Library does not contain cross section for " << name
<< " at or near " << std::round(temp_desired / K_BOLTZMANN) << "K.";
fatal_error(msg);
}
}
break;
case TEMPERATURE_INTERPOLATION:
// Get a list of bounding temperatures for each actual temperature
// present in the model
for (int k = 0; k < micros[m]->kTs.size() - 1; k++) {
for (int k = 0; k < micros[m]->kTs.shape()[0] - 1; k++) {
if ((micros[m]->kTs[k] <= temp_desired) &&
(temp_desired < micros[m]->kTs[k + 1])) {
micro_t[m] = k;
@ -394,8 +390,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
int num_interp_points = 2;
if (method == TEMPERATURE_NEAREST) num_interp_points = 1;
for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
double_1dvec interp(micros.size());
double_1dvec temp_indices(micros.size());
std::vector<double> interp(micros.size());
std::vector<double> temp_indices(micros.size());
for (int m = 0; m < micros.size(); m++) {
interp[m] = (1. - micro_t_interp[m]) * atom_densities[m];
temp_indices[m] = micro_t[m] + interp_point;
@ -409,8 +405,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
//==============================================================================
void
Mgxs::combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t)
Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scalars,
const std::vector<int>& micro_ts, int this_t)
{
// Build the vector of pointers to the xs objects within micros
std::vector<XsData*> those_xs(micros.size());
@ -441,19 +437,19 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs_t->total[a][gin];
val = xs_t->total(a, gin);
break;
case MG_GET_XS_NU_FISSION:
val = fissionable ? xs_t->nu_fission[a][gin] : 0.;
val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
break;
case MG_GET_XS_ABSORPTION:
val = xs_t->absorption[a][gin];
val = xs_t->absorption(a, gin);;
break;
case MG_GET_XS_FISSION:
val = fissionable ? xs_t->fission[a][gin] : 0.;
val = fissionable ? xs_t->fission(a, gin) : 0.;
break;
case MG_GET_XS_KAPPA_FISSION:
val = fissionable ? xs_t->kappa_fission[a][gin] : 0.;
val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
@ -462,16 +458,16 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_PROMPT_NU_FISSION:
val = fissionable ? xs_t->prompt_nu_fission[a][gin] : 0.;
val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
break;
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs_t->delayed_nu_fission[a][gin][*dg];
val = xs_t->delayed_nu_fission(a, *dg, gin);
} else {
val = 0.;
for (auto& num : xs_t->delayed_nu_fission[a][gin]) {
val += num;
for (int d = 0; d < xs_t->delayed_nu_fission.shape()[2]; d++) {
val += xs_t->delayed_nu_fission(a, d, gin);
}
}
} else {
@ -481,12 +477,12 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs_t->chi_prompt[a][gin][*gout];
val = xs_t->chi_prompt(a, gin, *gout);
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs_t->chi_prompt[a][gin]) {
val += num;
for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) {
val += xs_t->chi_prompt(a, gin, g);
}
}
} else {
@ -497,21 +493,21 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs_t->chi_delayed[a][gin][*gout][*dg];
val = xs_t->chi_delayed(a, *dg, gin, *gout);
} else {
val = xs_t->chi_delayed[a][gin][*gout][0];
val = xs_t->chi_delayed(a, 0, gin, *gout);
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
val += xs_t->chi_delayed[a][gin][i][*dg];
for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
val += xs_t->delayed_nu_fission(a, *dg, gin, g);
}
} else {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
for (auto& num : xs_t->chi_delayed[a][gin][i]) {
val += num;
for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) {
val += xs_t->delayed_nu_fission(a, d, gin, g);
}
}
}
@ -521,13 +517,13 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
}
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs_t->inverse_velocity[a][gin];
val = xs_t->inverse_velocity(a, gin);
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs_t->decay_rate[a][*dg + 1];
val = xs_t->decay_rate(a, *dg + 1);
} else {
val = xs_t->decay_rate[a][0];
val = xs_t->decay_rate(a, 0);
}
break;
default:
@ -548,11 +544,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
int tid = 0;
#endif
XsData* xs_t = &xs[cache[tid].t];
double nu_fission = xs_t->nu_fission[cache[tid].a][gin];
double nu_fission = xs_t->nu_fission(cache[tid].a, gin);
// Find the probability of having a prompt neutron
double prob_prompt =
xs_t->prompt_nu_fission[cache[tid].a][gin];
xs_t->prompt_nu_fission(cache[tid].a, gin);
// sample random numbers
double xi_pd = prn() * nu_fission;
@ -568,10 +564,10 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_prompt[cache[tid].a][gin][gout];
xs_t->chi_prompt(cache[tid].a, gin, gout);
while (prob_gout < xi_gout) {
gout++;
prob_gout += xs_t->chi_prompt[cache[tid].a][gin][gout];
prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout);
}
} else {
@ -582,7 +578,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs_t->delayed_nu_fission[cache[tid].a][gin][dg];
xs_t->delayed_nu_fission(cache[tid].a, dg, gin);
}
// adjust dg in case of round-off error
@ -591,11 +587,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
while (prob_gout < xi_gout) {
gout++;
prob_gout +=
xs_t->chi_delayed[cache[tid].a][gin][gout][dg];
xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
}
}
}
@ -630,10 +626,10 @@ Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3],
set_temperature_index(sqrtkT);
set_angle_index(uvw);
XsData* xs_t = &xs[cache[tid].t];
total_xs = xs_t->total[cache[tid].a][gin];
abs_xs = xs_t->absorption[cache[tid].a][gin];
total_xs = xs_t->total(cache[tid].a, gin);
abs_xs = xs_t->absorption(cache[tid].a, gin);
nu_fiss_xs = fissionable ? xs_t->nu_fission[cache[tid].a][gin] : 0.;
nu_fiss_xs = fissionable ? xs_t->nu_fission(cache[tid].a, gin) : 0.;
}
//==============================================================================
@ -662,17 +658,7 @@ Mgxs::set_temperature_index(double sqrtkT)
int tid = 0;
#endif
if (sqrtkT != cache[tid].sqrtkT) {
double kT = sqrtkT * sqrtkT;
// initialize vector for storage of the differences
std::valarray<double> temp_diff(kTs.data(), kTs.size());
// Find the minimum difference of kT and kTs
temp_diff = std::abs(temp_diff - kT);
cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// store this temperature as the last one used
cache[tid].t = xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
cache[tid].sqrtkT = sqrtkT;
}
}

6
src/mgxs_interface.cpp
View file

@ -19,7 +19,7 @@ add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
int& method)
{
// Convert temps to a vector for the from_hdf5 function
double_1dvec temperature(temps, temps + n_temps);
std::vector<double> temperature(temps, temps + n_temps);
write_message("Loading " + std::string(name) + " data...", 6);
@ -60,10 +60,10 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
{
if (n_temps > 0) {
// // Convert temps to a vector
double_1dvec temperature(temps, temps + n_temps);
std::vector<double> temperature(temps, temps + n_temps);
// Convert atom_densities to a vector
double_1dvec atom_densities_vec(atom_densities,
std::vector<double> atom_densities_vec(atom_densities,
atom_densities + n_nuclides);
// Build array of pointers to nuclides_MG's Mgxs objects needed for this

206
src/scattdata.cpp
View file

@ -4,6 +4,8 @@
#include <numeric>
#include <cmath>
#include "xtensor/xbuilder.hpp"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -16,11 +18,11 @@ namespace openmc {
//==============================================================================
void
ScattData::base_init(int order, const int_1dvec& in_gmin,
const int_1dvec& in_gmax, const double_2dvec& in_energy,
ScattData::base_init(int order, const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult)
{
int groups = in_energy.size();
size_t groups = in_energy.size();
gmin = in_gmin;
gmax = in_gmax;
@ -51,18 +53,17 @@ ScattData::base_init(int order, const int_1dvec& in_gmin,
//==============================================================================
void
ScattData::base_combine(int max_order,
const std::vector<ScattData*>& those_scatts, const double_1dvec& scalars,
int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult,
ScattData::base_combine(size_t max_order,
const std::vector<ScattData*>& those_scatts, const std::vector<double>& scalars,
xt::xtensor<int, 1>& in_gmin, xt::xtensor<int, 1>& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter)
{
int groups = those_scatts[0] -> energy.size();
size_t groups = those_scatts[0] -> energy.size();
// Now allocate and zero our storage spaces
double_3dvec this_matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(max_order, 0.)));
double_2dvec mult_numer(groups, double_1dvec(groups, 0.));
double_2dvec mult_denom(groups, double_1dvec(groups, 0.));
xt::xtensor<double, 3> this_matrix({groups, groups, max_order}, 0.);
xt::xtensor<double, 2> mult_numer({groups, groups}, 0.);
xt::xtensor<double, 2> mult_denom({groups, groups}, 0.);
// Build the dense scattering and multiplicity matrices
// Get the multiplicity_matrix
@ -80,26 +81,26 @@ ScattData::base_combine(int max_order,
ScattData* that = those_scatts[i];
// Build the dense matrix for that object
double_3dvec that_matrix = that->get_matrix(max_order);
xt::xtensor<double, 3> that_matrix = that->get_matrix(max_order);
// Now add that to this for the scattering and multiplicity
for (int gin = 0; gin < groups; gin++) {
// Only spend time adding that's gmin to gmax data since the rest will
// be zeros
int i_gout = 0;
for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) {
for (int gout = that->gmin(gin); gout <= that->gmax(gin); gout++) {
// Do the scattering matrix
for (int l = 0; l < max_order; l++) {
this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l];
this_matrix(gin, gout, l) += scalars[i] * that_matrix(gin, gout, l);
}
// Incorporate that's contribution to the multiplicity matrix data
double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout];
mult_numer[gin][gout] += scalars[i] * nuscatt;
double nuscatt = that->scattxs(gin) * that->energy[gin][i_gout];
mult_numer(gin, gout) += scalars[i] * nuscatt;
if (that->mult[gin][i_gout] > 0.) {
mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout];
mult_denom(gin, gout) += scalars[i] * nuscatt / that->mult[gin][i_gout];
} else {
mult_denom[gin][gout] += scalars[i];
mult_denom(gin, gout) += scalars[i];
}
i_gout++;
}
@ -107,16 +108,8 @@ ScattData::base_combine(int max_order,
}
// Combine mult_numer and mult_denom into the combined multiplicity matrix
double_2dvec this_mult(groups, double_1dvec(groups, 1.));
for (int gin = 0; gin < groups; gin++) {
for (int gout = 0; gout < groups; gout++) {
if (mult_denom[gin][gout] > 0.) {
this_mult[gin][gout] = mult_numer[gin][gout] / mult_denom[gin][gout];
}
}
}
mult_numer.clear();
mult_denom.clear();
xt::xtensor<double, 2> this_mult({groups, groups}, 1.);
this_mult = xt::nan_to_num(mult_numer / mult_denom);
// We have the data, now we need to convert to a jagged array and then use
// the initialize function to store it on the object.
@ -125,8 +118,8 @@ ScattData::base_combine(int max_order,
int gmin_;
for (gmin_ = 0; gmin_ < groups; gmin_++) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmin_].size(); l++) {
if (this_matrix[gin][gmin_][l] != 0.) {
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmin_, l) != 0.) {
non_zero = true;
break;
}
@ -136,8 +129,8 @@ ScattData::base_combine(int max_order,
int gmax_;
for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) {
bool non_zero = false;
for (int l = 0; l < this_matrix[gin][gmax_].size(); l++) {
if (this_matrix[gin][gmax_][l] != 0.) {
for (int l = 0; l < this_matrix.shape()[2]; l++) {
if (this_matrix(gin, gmax_, l) != 0.) {
non_zero = true;
break;
}
@ -160,8 +153,11 @@ ScattData::base_combine(int max_order,
sparse_mult[gin].resize(gmax_ - gmin_ + 1);
int i_gout = 0;
for (int gout = gmin_; gout <= gmax_; gout++) {
sparse_scatter[gin][i_gout] = this_matrix[gin][gout];
sparse_mult[gin][i_gout] = this_mult[gin][gout];
sparse_scatter[gin][i_gout].resize(this_matrix.shape()[2]);
for (int l = 0; l < this_matrix.shape()[2]; l++) {
sparse_scatter[gin][i_gout][l] = this_matrix(gin, gout, l);
}
sparse_mult[gin][i_gout] = this_mult(gin, gout);
i_gout++;
}
}
@ -241,21 +237,21 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
//==============================================================================
void
ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataLegendre::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
size_t groups = coeffs.size();
size_t order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Get the scattering cross section value by summing the un-normalized P0
// coefficient in the variable matrix over all outgoing groups.
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
int num_groups = in_gmax[gin] - in_gmin[gin] + 1;
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < num_groups; i_gout++) {
scattxs[gin] += matrix[gin][i_gout][0];
}
@ -301,7 +297,7 @@ ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
void
ScattDataLegendre::update_max_val()
{
int groups = max_val.size();
size_t groups = max_val.size();
// Step through the polynomial with fixed number of points to identify the
// maximal value
int Nmu = 1001;
@ -384,25 +380,25 @@ ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
void
ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
const std::vector<double>& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = 0;
size_t max_order = 0;
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataLegendre* that = dynamic_cast<ScattDataLegendre*>(those_scatts[i]);
if (!that) {
fatal_error("Cannot combine the ScattData objects!");
}
int that_order = that->get_order();
size_t that_order = that->get_order();
if (that_order > max_order) max_order = that_order;
}
max_order++; // Add one since this is a Legendre
int groups = those_scatts[0] -> energy.size();
size_t groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups}, 0);
xt::xtensor<int, 1> in_gmax({groups}, 0);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -418,20 +414,19 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
double_3dvec
ScattDataLegendre::get_matrix(int max_order)
xt::xtensor<double, 3>
ScattDataLegendre::get_matrix(size_t max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
int order_dim = max_order + 1;
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
size_t groups = energy.size();
size_t order_dim = max_order + 1;
xt::xtensor<double, 3> matrix({groups, groups, order_dim}, 0.);
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
dist[gin][i_gout][l];
}
}
@ -444,20 +439,20 @@ ScattDataLegendre::get_matrix(int max_order)
//==============================================================================
void
ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataHistogram::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
size_t groups = coeffs.size();
size_t order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Get the scattering cross section value by summing the distribution
// over all the histogram bins in angle and outgoing energy groups
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
scattxs[gin] += std::accumulate(matrix[gin][i_gout].begin(),
matrix[gin][i_gout].end(), 0.);
@ -484,12 +479,8 @@ ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Build the angular distribution mu values
mu = double_1dvec(order);
mu = xt::linspace(-1., 1., order + 1);
dmu = 2. / order;
mu[0] = -1.;
for (int imu = 1; imu < order; imu++) {
mu[imu] = -1. + imu * dmu;
}
// Calculate f(mu) and integrate it so we can avoid rejection sampling
fmu.resize(groups);
@ -534,7 +525,7 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu)
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
imu = this->mu.shape()[0] - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
@ -560,7 +551,6 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
if (xi < dist[gin][i_gout][0]) {
imu = 0;
} else {
// TODO lower_bound? + 1?
imu = std::upper_bound(dist[gin][i_gout].begin(),
dist[gin][i_gout].end(), xi) -
dist[gin][i_gout].begin();
@ -581,21 +571,20 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataHistogram::get_matrix(int max_order)
xt::xtensor<double, 3>
ScattDataHistogram::get_matrix(size_t max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
size_t groups = energy.size();
// We ignore the requested order for Histogram and Tabular representations
int order_dim = get_order();
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
size_t order_dim = get_order();
xt::xtensor<double, 3> matrix({groups, groups, order_dim}, 0);
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
@ -607,10 +596,10 @@ ScattDataHistogram::get_matrix(int max_order)
void
ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
const std::vector<double>& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = those_scatts[0]->get_order();
size_t max_order = those_scatts[0]->get_order();
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataHistogram* that = dynamic_cast<ScattDataHistogram*>(those_scatts[i]);
@ -622,10 +611,10 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
}
}
int groups = those_scatts[0] -> energy.size();
size_t groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups}, 0);
xt::xtensor<int, 1> in_gmax({groups}, 0);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -633,7 +622,7 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -644,29 +633,24 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
void
ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
ScattDataTabular::init(const xt::xtensor<int, 1>& in_gmin,
const xt::xtensor<int, 1>& in_gmax, const double_2dvec& in_mult,
const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
size_t groups = coeffs.size();
size_t order = coeffs[0][0].size();
// make a copy of coeffs that we can use to both extract data and normalize
double_3dvec matrix = coeffs;
// Build the angular distribution mu values
mu = double_1dvec(order);
mu = xt::linspace(-1., 1., order);
dmu = 2. / (order - 1);
mu[0] = -1.;
for (int imu = 1; imu < order - 1; imu++) {
mu[imu] = -1. + imu * dmu;
}
mu[order - 1] = 1.;
// Get the scattering cross section value by integrating the distribution
// over all mu points and then combining over all outgoing groups
scattxs.resize(groups);
scattxs = xt::zeros<double>({groups});
for (int gin = 0; gin < groups; gin++) {
scattxs[gin] = 0.;
for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) {
for (int imu = 1; imu < order; imu++) {
scattxs[gin] += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] +
@ -743,7 +727,7 @@ ScattDataTabular::calc_f(int gin, int gout, double mu)
int imu;
if (mu == 1.) {
// use size -2 to have the index one before the end
imu = this->mu.size() - 2;
imu = this->mu.shape()[0] - 2;
} else {
imu = std::floor((mu + 1.) / dmu + 1.) - 1;
}
@ -764,7 +748,7 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
sample_energy(gin, gout, i_gout);
// Determine the outgoing cosine bin
int NP = this->mu.size();
int NP = this->mu.shape()[0];
double xi = prn();
double c_k = dist[gin][i_gout][0];
@ -804,21 +788,20 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataTabular::get_matrix(int max_order)
xt::xtensor<double, 3>
ScattDataTabular::get_matrix(size_t max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
size_t groups = energy.size();
// We ignore the requested order for Histogram and Tabular representations
int order_dim = get_order();
double_3dvec matrix = double_3dvec(groups, double_2dvec(groups,
double_1dvec(order_dim, 0.)));
size_t order_dim = get_order();
xt::xtensor<double, 3> matrix({groups, groups, order_dim}, 0.);
for (int gin = 0; gin < groups; gin++) {
for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) {
int gout = i_gout + gmin[gin];
for (int l = 0; l < order_dim; l++) {
matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] *
matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] *
fmu[gin][i_gout][l];
}
}
@ -830,10 +813,10 @@ ScattDataTabular::get_matrix(int max_order)
void
ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars)
const std::vector<double>& scalars)
{
// Find the max order in the data set and make sure we can combine the sets
int max_order = those_scatts[0]->get_order();
size_t max_order = those_scatts[0]->get_order();
for (int i = 0; i < those_scatts.size(); i++) {
// Lets also make sure these items are combineable
ScattDataTabular* that = dynamic_cast<ScattDataTabular*>(those_scatts[i]);
@ -845,10 +828,10 @@ ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
}
}
int groups = those_scatts[0] -> energy.size();
size_t groups = those_scatts[0] -> energy.size();
int_1dvec in_gmin(groups);
int_1dvec in_gmax(groups);
xt::xtensor<int, 1> in_gmin({groups}, 0);
xt::xtensor<int, 1> in_gmax({groups}, 0);
double_3dvec sparse_scatter(groups);
double_2dvec sparse_mult(groups);
@ -856,7 +839,7 @@ ScattDataTabular::combine(const std::vector<ScattData*>& those_scatts,
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -885,16 +868,11 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
tab.scattxs = leg.scattxs;
// Build mu and dmu
tab.mu = double_1dvec(n_mu);
tab.mu = xt::linspace(-1., 1., n_mu);
tab.dmu = 2. / (n_mu - 1);
tab.mu[0] = -1.;
for (int imu = 1; imu < n_mu - 1; imu++) {
tab.mu[imu] = -1. + imu * tab.dmu;
}
tab.mu[n_mu - 1] = 1.;
// Calculate f(mu) and integrate it so we can avoid rejection sampling
int groups = tab.energy.size();
size_t groups = tab.energy.size();
tab.fmu.resize(groups);
for (int gin = 0; gin < groups; gin++) {
int num_groups = tab.gmax[gin] - tab.gmin[gin] + 1;

864
src/xsdata.cpp
View file

File diff suppressed because it is too large Load diff

BIN
tests/1d_mgxs.h5
View file

Binary file not shown.

103
tests/regression_tests/mg_basic/inputs_true.dat
View file

@ -1,97 +1,64 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="4" material="4" region="1 -2 3 -4 8 -9" universe="1" />
<cell id="5" material="5" region="1 -2 3 -4 9 -10" universe="1" />
<cell id="6" material="6" region="1 -2 3 -4 10 -11" universe="1" />
<cell id="7" material="7" region="1 -2 3 -4 11 -12" universe="1" />
<cell id="8" material="8" region="1 -2 3 -4 12 -13" universe="1" />
<cell id="9" material="9" region="1 -2 3 -4 13 -14" universe="1" />
<cell id="10" material="10" region="1 -2 3 -4 14 -15" universe="1" />
<cell id="11" material="11" region="1 -2 3 -4 15 -16" universe="1" />
<cell id="12" material="12" region="1 -2 3 -4 16 -17" universe="1" />
<cell id="1" material="1" region="1 -2" universe="0" />
<cell id="2" material="2" region="2 -3" universe="0" />
<cell id="3" material="3" region="3 -4" universe="0" />
<cell id="4" material="4" region="4 -5" universe="0" />
<cell id="5" material="5" region="5 -6" universe="0" />
<cell id="6" material="6" region="6 -7" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="0.4167" id="6" type="z-plane" />
<surface coeffs="0.8334" id="7" type="z-plane" />
<surface coeffs="1.2501" id="8" type="z-plane" />
<surface coeffs="1.6668" id="9" type="z-plane" />
<surface coeffs="2.0835" id="10" type="z-plane" />
<surface coeffs="2.5002" id="11" type="z-plane" />
<surface coeffs="2.9169" id="12" type="z-plane" />
<surface coeffs="3.3336" id="13" type="z-plane" />
<surface coeffs="3.7503" id="14" type="z-plane" />
<surface coeffs="4.167" id="15" type="z-plane" />
<surface coeffs="4.5837" id="16" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="17" type="z-plane" />
<surface coeffs="154.90833333333333" id="2" type="x-plane" />
<surface coeffs="309.81666666666666" id="3" type="x-plane" />
<surface coeffs="464.725" id="4" type="x-plane" />
<surface coeffs="619.6333333333333" id="5" type="x-plane" />
<surface coeffs="774.5416666666666" id="6" type="x-plane" />
<surface boundary="vacuum" coeffs="929.45" id="7" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<cross_sections>2g.h5</cross_sections>
<material id="1" name="base leg">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang" />
<macroscopic name="mat_1" />
</material>
<material id="2" name="2">
<material id="2" name="base tab">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang_mu" />
<macroscopic name="mat_2" />
</material>
<material id="3" name="3">
<material id="3" name="base hist">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
<macroscopic name="mat_3" />
</material>
<material id="4" name="4">
<material id="4" name="base matrix">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso_mu" />
<macroscopic name="mat_4" />
</material>
<material id="5" name="5">
<material id="5" name="base ang">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang" />
<macroscopic name="mat_5" />
</material>
<material id="6" name="6">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang_mu" />
</material>
<material id="7" name="7">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="8" name="8">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso_mu" />
</material>
<material id="9" name="9">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang" />
</material>
<material id="10" name="10">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang_mu" />
</material>
<material id="11" name="11">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
</material>
<material id="12" name="12">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso_mu" />
<material id="6" name="micro">
<density units="sum" />
<nuclide ao="0.5" name="mat_1" />
<nuclide ao="0.5" name="mat_6" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
<parameters>0.0 -1000.0 -1000.0 154.90833333333333 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>

2
tests/regression_tests/mg_basic/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.073147E+00 1.602384E-02
1.005345E+00 1.109180E-02

90
tests/regression_tests/mg_basic/test.py
View file

@ -1,9 +1,97 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Make a version of mat-1 which has a tabular representation of the
# scattering vice Legendre with 33 points
mat_2 = mat_1.convert_scatter_format('tabular', 33)
mat_2.name = 'mat_2'
mg_cross_sections_file.add_xsdata(mat_2)
# Make a version of mat-1 which has a histogram representation of the
# scattering vice Legendre with 33 bins
mat_3 = mat_1.convert_scatter_format('histogram', 33)
mat_3.name = 'mat_3'
mg_cross_sections_file.add_xsdata(mat_3)
# Make a version which uses a fission matrix vice chi & nu-fission
mat_4 = openmc.XSdata('mat_4', groups)
mat_4.order = 1
mat_4.set_nu_fission(np.outer(np.multiply(nu, fiss), chi))
mat_4.set_absorption(absorption)
mat_4.set_scatter_matrix(scatter)
mat_4.set_total(total)
mg_cross_sections_file.add_xsdata(mat_4)
# Make an angle-dependent version of mat_1 with 2 polar and 2 azim. angles
mat_5 = mat_1.convert_representation('angle', 2, 2)
mat_5.name = 'mat_5'
mg_cross_sections_file.add_xsdata(mat_5)
# Make a copy of mat_1 for testing microscopic cross sections
mat_6 = openmc.XSdata('mat_6', groups)
mat_6.order = 1
mat_6.set_nu_fission(np.multiply(nu, fiss))
mat_6.set_absorption(absorption)
mat_6.set_scatter_matrix(scatter)
mat_6.set_total(total)
mat_6.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_6)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_basic():
model = slab_mg()
create_library()
mat_names = ['base leg', 'base tab', 'base hist', 'base matrix',
'base ang', 'micro']
model = slab_mg(num_regions=6, mat_names=mat_names)
# Modify the last material to be a microscopic combination of nuclides
model.materials[-1] = openmc.Material(name='micro', material_id=6)
model.materials[-1].set_density("sum")
model.materials[-1].add_nuclide("mat_1", 0.5)
model.materials[-1].add_nuclide("mat_6", 0.5)
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

0
tests/regression_tests/mg_nuclide/__init__.py → tests/regression_tests/mg_basic_delayed/__init__.py
View file

63
tests/regression_tests/mg_basic_delayed/inputs_true.dat Normal file
View file

@ -0,0 +1,63 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2" universe="0" />
<cell id="2" material="2" region="2 -3" universe="0" />
<cell id="3" material="3" region="3 -4" universe="0" />
<cell id="4" material="4" region="4 -5" universe="0" />
<cell id="5" material="5" region="5 -6" universe="0" />
<cell id="6" material="6" region="6 -7" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface coeffs="154.90833333333333" id="2" type="x-plane" />
<surface coeffs="309.81666666666666" id="3" type="x-plane" />
<surface coeffs="464.725" id="4" type="x-plane" />
<surface coeffs="619.6333333333333" id="5" type="x-plane" />
<surface coeffs="774.5416666666666" id="6" type="x-plane" />
<surface boundary="vacuum" coeffs="929.45" id="7" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>2g.h5</cross_sections>
<material id="1" name="vec beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_1" />
</material>
<material id="2" name="vec no beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_2" />
</material>
<material id="3" name="matrix beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_3" />
</material>
<material id="4" name="matrix no beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_4" />
</material>
<material id="5" name="vec group beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_5" />
</material>
<material id="6" name="matrix group beta">
<density units="macro" value="1.0" />
<macroscopic name="mat_6" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 -1000.0 -1000.0 154.90833333333333 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>

2
tests/regression_tests/mg_basic_delayed/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.003463E+00 2.173155E-02

131
tests/regression_tests/mg_basic_delayed/test.py Normal file
View file

@ -0,0 +1,131 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
n_dg = 2
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.num_delayed_groups = n_dg
beta = np.array([0.003, 0.003])
one_m_beta = 1. - np.sum(beta)
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
# Make the base data that uses chi & nu-fission vectors with a beta
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.num_delayed_groups = 2
mat_1.set_beta(beta)
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Make a version that uses prompt and delayed version of nufiss and chi
mat_2 = openmc.XSdata('mat_2', groups)
mat_2.order = 1
mat_2.num_delayed_groups = 2
mat_2.set_prompt_nu_fission(one_m_beta * np.multiply(nu, fiss))
delay_nu_fiss = np.zeros((n_dg, groups.num_groups))
for dg in range(n_dg):
for g in range(groups.num_groups):
delay_nu_fiss[dg, g] = beta[dg] * nu[g] * fiss[g]
mat_2.set_delayed_nu_fission(delay_nu_fiss)
mat_2.set_absorption(absorption)
mat_2.set_scatter_matrix(scatter)
mat_2.set_total(total)
mat_2.set_chi_prompt(chi)
mat_2.set_chi_delayed(np.stack([chi] * n_dg))
mg_cross_sections_file.add_xsdata(mat_2)
# Make a version that uses a nu-fission matrix with a beta
mat_3 = openmc.XSdata('mat_3', groups)
mat_3.order = 1
mat_3.num_delayed_groups = 2
mat_3.set_beta(beta)
mat_3.set_nu_fission(np.outer(np.multiply(nu, fiss), chi))
mat_3.set_absorption(absorption)
mat_3.set_scatter_matrix(scatter)
mat_3.set_total(total)
mg_cross_sections_file.add_xsdata(mat_3)
# Make a version that uses prompt and delayed version of the nufiss matrix
mat_4 = openmc.XSdata('mat_4', groups)
mat_4.order = 1
mat_4.num_delayed_groups = 2
mat_4.set_prompt_nu_fission(one_m_beta *
np.outer(np.multiply(nu, fiss), chi))
delay_nu_fiss = np.zeros((n_dg, groups.num_groups, groups.num_groups))
for dg in range(n_dg):
for g in range(groups.num_groups):
for go in range(groups.num_groups):
delay_nu_fiss[dg, g, go] = beta[dg] * nu[g] * fiss[g] * chi[go]
mat_4.set_delayed_nu_fission(delay_nu_fiss)
mat_4.set_absorption(absorption)
mat_4.set_scatter_matrix(scatter)
mat_4.set_total(total)
mg_cross_sections_file.add_xsdata(mat_4)
# Make the base data that uses chi & nu-fiss vectors with a group-wise beta
mat_5 = openmc.XSdata('mat_5', groups)
mat_5.order = 1
mat_5.num_delayed_groups = 2
mat_5.set_beta(np.stack([beta] * groups.num_groups))
mat_5.set_nu_fission(np.multiply(nu, fiss))
mat_5.set_absorption(absorption)
mat_5.set_scatter_matrix(scatter)
mat_5.set_total(total)
mat_5.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_5)
# Make a version that uses a nu-fission matrix with a group-wise beta
mat_6 = openmc.XSdata('mat_6', groups)
mat_6.order = 1
mat_6.num_delayed_groups = 2
mat_6.set_beta(np.stack([beta] * groups.num_groups))
mat_6.set_nu_fission(np.outer(np.multiply(nu, fiss), chi))
mat_6.set_absorption(absorption)
mat_6.set_scatter_matrix(scatter)
mat_6.set_total(total)
mg_cross_sections_file.add_xsdata(mat_6)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_basic_delayed():
create_library()
model = slab_mg(num_regions=6, mat_names=['vec beta', 'vec no beta',
'matrix beta', 'matrix no beta',
'vec group beta',
'matrix group beta'])
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

2
tests/regression_tests/mg_convert/test.py
View file

@ -17,7 +17,7 @@ def build_mgxs_library(convert):
# Instantiate the energy group data
groups = openmc.mgxs.EnergyGroups(group_edges=[1e-5, 0.625, 20.0e6])
# Instantiate the 7-group (C5G7) cross section data
# Instantiate the 2-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2', groups)
uo2_xsdata.order = 2
uo2_xsdata.set_total([2., 2.])

33
tests/regression_tests/mg_legendre/inputs_true.dat
View file

@ -1,44 +1,31 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="1" material="1" region="1 -2" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="1.6667" id="6" type="z-plane" />
<surface coeffs="3.3334" id="7" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="8" type="z-plane" />
<surface boundary="vacuum" coeffs="929.45" id="2" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<cross_sections>2g.h5</cross_sections>
<material id="1" name="mat_1">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
</material>
<material id="2" name="2">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="3" name="3">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
<macroscopic name="mat_1" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
<parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>

2
tests/regression_tests/mg_legendre/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.110122E+00 2.549637E-02
9.934975E-01 2.679669E-02

46
tests/regression_tests/mg_legendre/test.py
View file

@ -1,10 +1,54 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_legendre():
model = slab_mg(reps=['iso'])
create_library()
model = slab_mg()
model.settings.tabular_legendre = {'enable': False}
harness = PyAPITestHarness('statepoint.10.h5', model)

36
tests/regression_tests/mg_max_order/inputs_true.dat
View file

@ -1,44 +1,34 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="1" material="1" region="1 -2" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="1.6667" id="6" type="z-plane" />
<surface coeffs="3.3334" id="7" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="8" type="z-plane" />
<surface boundary="vacuum" coeffs="929.45" id="2" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<cross_sections>2g.h5</cross_sections>
<material id="1" name="mat_1">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
</material>
<material id="2" name="2">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="3" name="3">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
<macroscopic name="mat_1" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
<parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<max_order>1</max_order>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>

2
tests/regression_tests/mg_max_order/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.074551E+00 1.871525E-02
9.934975E-01 2.679669E-02

47
tests/regression_tests/mg_max_order/test.py
View file

@ -1,10 +1,55 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694, 0.003], [0.004555, -0.0003972, 0.00002]],
[[0.00000, 0.00000, 0.000], [0.424100, 0.05439000, 0.0025]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 2
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_max_order():
model = slab_mg(reps=['iso'])
create_library()
model = slab_mg()
model.settings.max_order = 1
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

97
tests/regression_tests/mg_nuclide/inputs_true.dat
View file

@ -1,97 +0,0 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="4" material="4" region="1 -2 3 -4 8 -9" universe="1" />
<cell id="5" material="5" region="1 -2 3 -4 9 -10" universe="1" />
<cell id="6" material="6" region="1 -2 3 -4 10 -11" universe="1" />
<cell id="7" material="7" region="1 -2 3 -4 11 -12" universe="1" />
<cell id="8" material="8" region="1 -2 3 -4 12 -13" universe="1" />
<cell id="9" material="9" region="1 -2 3 -4 13 -14" universe="1" />
<cell id="10" material="10" region="1 -2 3 -4 14 -15" universe="1" />
<cell id="11" material="11" region="1 -2 3 -4 15 -16" universe="1" />
<cell id="12" material="12" region="1 -2 3 -4 16 -17" universe="1" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="0.4167" id="6" type="z-plane" />
<surface coeffs="0.8334" id="7" type="z-plane" />
<surface coeffs="1.2501" id="8" type="z-plane" />
<surface coeffs="1.6668" id="9" type="z-plane" />
<surface coeffs="2.0835" id="10" type="z-plane" />
<surface coeffs="2.5002" id="11" type="z-plane" />
<surface coeffs="2.9169" id="12" type="z-plane" />
<surface coeffs="3.3336" id="13" type="z-plane" />
<surface coeffs="3.7503" id="14" type="z-plane" />
<surface coeffs="4.167" id="15" type="z-plane" />
<surface coeffs="4.5837" id="16" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="17" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_ang" />
</material>
<material id="2" name="2">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_ang_mu" />
</material>
<material id="3" name="3">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_iso" />
</material>
<material id="4" name="4">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_iso_mu" />
</material>
<material id="5" name="5">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_ang" />
</material>
<material id="6" name="6">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_ang_mu" />
</material>
<material id="7" name="7">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_iso" />
</material>
<material id="8" name="8">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_iso_mu" />
</material>
<material id="9" name="9">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_ang" />
</material>
<material id="10" name="10">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_ang_mu" />
</material>
<material id="11" name="11">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_iso" />
</material>
<material id="12" name="12">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_iso_mu" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

2
tests/regression_tests/mg_nuclide/results_true.dat
View file

@ -1,2 +0,0 @@
k-combined:
1.073147E+00 1.602384E-02

9
tests/regression_tests/mg_nuclide/test.py
View file

@ -1,9 +0,0 @@
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def test_mg_nuclide():
model = slab_mg(as_macro=False)
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

90
tests/regression_tests/mg_survival_biasing/inputs_true.dat
View file

@ -1,98 +1,34 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="4" material="4" region="1 -2 3 -4 8 -9" universe="1" />
<cell id="5" material="5" region="1 -2 3 -4 9 -10" universe="1" />
<cell id="6" material="6" region="1 -2 3 -4 10 -11" universe="1" />
<cell id="7" material="7" region="1 -2 3 -4 11 -12" universe="1" />
<cell id="8" material="8" region="1 -2 3 -4 12 -13" universe="1" />
<cell id="9" material="9" region="1 -2 3 -4 13 -14" universe="1" />
<cell id="10" material="10" region="1 -2 3 -4 14 -15" universe="1" />
<cell id="11" material="11" region="1 -2 3 -4 15 -16" universe="1" />
<cell id="12" material="12" region="1 -2 3 -4 16 -17" universe="1" />
<cell id="1" material="1" region="1 -2" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="0.4167" id="6" type="z-plane" />
<surface coeffs="0.8334" id="7" type="z-plane" />
<surface coeffs="1.2501" id="8" type="z-plane" />
<surface coeffs="1.6668" id="9" type="z-plane" />
<surface coeffs="2.0835" id="10" type="z-plane" />
<surface coeffs="2.5002" id="11" type="z-plane" />
<surface coeffs="2.9169" id="12" type="z-plane" />
<surface coeffs="3.3336" id="13" type="z-plane" />
<surface coeffs="3.7503" id="14" type="z-plane" />
<surface coeffs="4.167" id="15" type="z-plane" />
<surface coeffs="4.5837" id="16" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="17" type="z-plane" />
<surface boundary="vacuum" coeffs="929.45" id="2" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<cross_sections>2g.h5</cross_sections>
<material id="1" name="mat_1">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang" />
</material>
<material id="2" name="2">
<density units="macro" value="1.0" />
<macroscopic name="uo2_ang_mu" />
</material>
<material id="3" name="3">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
</material>
<material id="4" name="4">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso_mu" />
</material>
<material id="5" name="5">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang" />
</material>
<material id="6" name="6">
<density units="macro" value="1.0" />
<macroscopic name="clad_ang_mu" />
</material>
<material id="7" name="7">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="8" name="8">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso_mu" />
</material>
<material id="9" name="9">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang" />
</material>
<material id="10" name="10">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_ang_mu" />
</material>
<material id="11" name="11">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
</material>
<material id="12" name="12">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso_mu" />
<macroscopic name="mat_1" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
<parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<survival_biasing>true</survival_biasing>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>

2
tests/regression_tests/mg_survival_biasing/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.080832E+00 1.336780E-02
9.979905E-01 6.207495E-03

45
tests/regression_tests/mg_survival_biasing/test.py
View file

@ -1,10 +1,55 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_survival_biasing():
create_library()
model = slab_mg()
model.settings.survival_biasing = True
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

122
tests/regression_tests/mg_tallies/inputs_true.dat
View file

@ -1,112 +1,48 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2 3 -4 5 -6" universe="1" />
<cell id="2" material="2" region="1 -2 3 -4 6 -7" universe="1" />
<cell id="3" material="3" region="1 -2 3 -4 7 -8" universe="1" />
<cell id="4" material="4" region="1 -2 3 -4 8 -9" universe="1" />
<cell id="5" material="5" region="1 -2 3 -4 9 -10" universe="1" />
<cell id="6" material="6" region="1 -2 3 -4 10 -11" universe="1" />
<cell id="7" material="7" region="1 -2 3 -4 11 -12" universe="1" />
<cell id="8" material="8" region="1 -2 3 -4 12 -13" universe="1" />
<cell id="9" material="9" region="1 -2 3 -4 13 -14" universe="1" />
<cell id="10" material="10" region="1 -2 3 -4 14 -15" universe="1" />
<cell id="11" material="11" region="1 -2 3 -4 15 -16" universe="1" />
<cell id="12" material="12" region="1 -2 3 -4 16 -17" universe="1" />
<cell id="1" material="1" region="1 -2" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="2" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="3" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="4" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="5" type="z-plane" />
<surface coeffs="0.4167" id="6" type="z-plane" />
<surface coeffs="0.8334" id="7" type="z-plane" />
<surface coeffs="1.2501" id="8" type="z-plane" />
<surface coeffs="1.6668" id="9" type="z-plane" />
<surface coeffs="2.0835" id="10" type="z-plane" />
<surface coeffs="2.5002" id="11" type="z-plane" />
<surface coeffs="2.9169" id="12" type="z-plane" />
<surface coeffs="3.3336" id="13" type="z-plane" />
<surface coeffs="3.7503" id="14" type="z-plane" />
<surface coeffs="4.167" id="15" type="z-plane" />
<surface coeffs="4.5837" id="16" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="17" type="z-plane" />
<surface boundary="vacuum" coeffs="929.45" id="2" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../../1d_mgxs.h5</cross_sections>
<material id="1" name="1">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_ang" />
</material>
<material id="2" name="2">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_ang_mu" />
</material>
<material id="3" name="3">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_iso" />
</material>
<material id="4" name="4">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="uo2_iso_mu" />
</material>
<material id="5" name="5">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_ang" />
</material>
<material id="6" name="6">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_ang_mu" />
</material>
<material id="7" name="7">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_iso" />
</material>
<material id="8" name="8">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="clad_iso_mu" />
</material>
<material id="9" name="9">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_ang" />
</material>
<material id="10" name="10">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_ang_mu" />
</material>
<material id="11" name="11">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_iso" />
</material>
<material id="12" name="12">
<density units="atom/b-cm" value="1.0" />
<nuclide ao="1.0" name="lwtr_iso_mu" />
<cross_sections>2g.h5</cross_sections>
<material id="1" name="mat_1">
<density units="macro" value="1.0" />
<macroscopic name="mat_1" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
<parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
</space>
</source>
<output>
<summary>false</summary>
</output>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<dimension>1 1 10</dimension>
<dimension>10 1 1</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>10 10 5</upper_right>
<upper_right>929.45 1000 1000</upper_right>
</mesh>
<filter id="5" type="mesh">
<bins>1</bins>
</filter>
<filter id="6" type="material">
<bins>1 2 3 4 5 6 7 8 9 10 11 12</bins>
<bins>1</bins>
</filter>
<filter id="1" type="energy">
<bins>0.0 20000000.0</bins>
@ -115,10 +51,10 @@
<bins>0.0 20000000.0</bins>
</filter>
<filter id="3" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
<bins>0.0 0.625 20000000.0</bins>
</filter>
<filter id="4" type="energyout">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
<bins>0.0 0.625 20000000.0</bins>
</filter>
<tally id="1">
<filters>5</filters>
@ -170,60 +106,60 @@
</tally>
<tally id="11">
<filters>5</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>analog</estimator>
</tally>
<tally id="12">
<filters>5</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="13">
<filters>6 1</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission scatter nu-scatter</scores>
<estimator>analog</estimator>
</tally>
<tally id="14">
<filters>6 1</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>collision</estimator>
</tally>
<tally id="15">
<filters>6 1</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="16">
<filters>6 1 2</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>scatter nu-scatter nu-fission</scores>
</tally>
<tally id="17">
<filters>6 3</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission scatter nu-scatter</scores>
<estimator>analog</estimator>
</tally>
<tally id="18">
<filters>6 3</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>collision</estimator>
</tally>
<tally id="19">
<filters>6 3</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>total absorption fission nu-fission</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="20">
<filters>6 3 4</filters>
<nuclides>uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu</nuclides>
<nuclides>mat_1</nuclides>
<scores>scatter nu-scatter nu-fission</scores>
</tally>
</tallies>

2
tests/regression_tests/mg_tallies/results_true.dat
View file

@ -1 +1 @@
9183f8b191f2e62334f992acd865d29e3f4e3f871a6df498e280fc4e2d91f2d2d20c732fbd75fa88e2e8c576f86e744f7655af6bb9da66e9b28b1009c8742899
41ea1f6b17c58a8141921af2f1d044eda93f3a9bca9463ee023af2e9865da613ace90fc8a25b42edde128ed827182ea9df0fe09d9b7887282d0ec092692cf717

51
tests/regression_tests/mg_tallies/test.py
View file

@ -1,23 +1,66 @@
import os
import numpy as np
import openmc
from openmc.examples import slab_mg
from tests.testing_harness import HashedPyAPITestHarness
def create_library():
# Instantiate the energy group data and file object
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
mg_cross_sections_file = openmc.MGXSLibrary(groups)
# Make the base, isotropic data
nu = [2.50, 2.50]
fiss = np.array([0.002817, 0.097])
capture = [0.008708, 0.02518]
absorption = np.add(capture, fiss)
scatter = np.array(
[[[0.31980, 0.06694], [0.004555, -0.0003972]],
[[0.00000, 0.00000], [0.424100, 0.05439000]]])
total = [0.33588, 0.54628]
chi = [1., 0.]
mat_1 = openmc.XSdata('mat_1', groups)
mat_1.order = 1
mat_1.set_nu_fission(np.multiply(nu, fiss))
mat_1.set_absorption(absorption)
mat_1.set_scatter_matrix(scatter)
mat_1.set_total(total)
mat_1.set_chi(chi)
mg_cross_sections_file.add_xsdata(mat_1)
# Write the file
mg_cross_sections_file.export_to_hdf5('2g.h5')
class MGXSTestHarness(HashedPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = '2g.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_tallies():
model = slab_mg(as_macro=False)
create_library()
model = slab_mg()
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh.type = 'regular'
mesh.dimension = [1, 1, 10]
mesh.dimension = [10, 1, 1]
mesh.lower_left = [0.0, 0.0, 0.0]
mesh.upper_right = [10, 10, 5]
mesh.upper_right = [929.45, 1000, 1000]
# Instantiate some tally filters
energy_filter = openmc.EnergyFilter([0.0, 20.0e6])
energyout_filter = openmc.EnergyoutFilter([0.0, 20.0e6])
energies = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
energies = [0.0, 0.625, 20.0e6]
matching_energy_filter = openmc.EnergyFilter(energies)
matching_eout_filter = openmc.EnergyoutFilter(energies)
mesh_filter = openmc.MeshFilter(mesh)