Multigroup Per-Thread Cache Conversion (#2591)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
John Tramm 2023-07-21 02:53:32 -05:00 committed by GitHub
parent 8101328d96
commit b5001ee12a
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 487 additions and 225 deletions

View file

@ -16,20 +16,6 @@
namespace openmc {
//==============================================================================
// Cache contains the cached data for an MGXS object
//==============================================================================
struct CacheData {
double sqrtkT; // last temperature corresponding to t
int t; // temperature index
int a; // angle index
// last angle that corresponds to a
double u;
double v;
double w;
};
//==============================================================================
// MGXS contains the mgxs data for a nuclide/material
//==============================================================================
@ -96,10 +82,9 @@ private:
bool equiv(const Mgxs& that);
public:
std::string name; // name of dataset, e.g., UO2
double awr; // atomic weight ratio
bool fissionable; // Is this fissionable
vector<CacheData> cache; // index and data cache
std::string name; // name of dataset, e.g., UO2
double awr; // atomic weight ratio
bool fissionable; // Is this fissionable
Mgxs() = default;
@ -135,13 +120,15 @@ public:
//! @param mu Cosine of the change-in-angle, for scattering quantities;
//! use nullptr if irrelevant.
//! @param dg delayed group index; use nullptr if irrelevant.
//! @param t Temperature index.
//! @param a Angle index.
//! @return Requested cross section value.
double get_xs(
MgxsType xstype, int gin, const int* gout, const double* mu, const int* dg);
double get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
const int* dg, int t, int a);
inline double get_xs(MgxsType xstype, int gin)
inline double get_xs(MgxsType xstype, int gin, int t, int a)
{
return get_xs(xstype, gin, nullptr, nullptr, nullptr);
return get_xs(xstype, gin, nullptr, nullptr, nullptr, t, a);
}
//! \brief Samples the fission neutron energy and if prompt or delayed.
@ -150,7 +137,10 @@ public:
//! @param dg Sampled delayed group index.
//! @param gout Sampled outgoing energy group.
//! @param seed Pseudorandom seed pointer
void sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed);
//! @param t Temperature index.
//! @param a Angle index.
void sample_fission_energy(
int gin, int& dg, int& gout, uint64_t* seed, int t, int a);
//! \brief Samples the outgoing energy and angle from a scatter event.
//!
@ -159,23 +149,37 @@ public:
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
//! @param seed Pseudorandom seed pointer.
//! @param t Temperature index.
//! @param a Angle index.
void sample_scatter(
int gin, int& gout, double& mu, double& wgt, uint64_t* seed);
int gin, int& gout, double& mu, double& wgt, uint64_t* seed, int t, int a);
//! \brief Calculates cross section quantities needed for tracking.
//!
//! @param p The particle whose attributes set which MGXS to get.
void calculate_xs(Particle& p);
//! \brief Sets the temperature index in cache given a temperature
//! \brief Sets the temperature index in the particle's cache.
//!
//! @param p Particle.
void set_temperature_index(Particle& p);
//! \brief Gets the temperature index given a temperature.
//!
//! @param sqrtkT Temperature of the material.
void set_temperature_index(double sqrtkT);
//! @return The temperature index corresponding to sqrtkT.
int get_temperature_index(double sqrtkT) const;
//! \brief Sets the angle index in cache given a direction
//! \brief Sets the angle index in the particle's cache.
//!
//! @param p Particle.
void set_angle_index(Particle& p);
//! \brief Gets the angle index given a direction.
//!
//! @param u Incoming particle direction.
void set_angle_index(Direction u);
//! @return The angle index corresponding to u.
int get_angle_index(const Direction& u) const;
//! \brief Provide const access to list of XsData held by this
const vector<XsData>& get_xsdata() const { return xs; }

View file

@ -170,6 +170,18 @@ struct MacroXS {
double pair_production; //!< macroscopic pair production xs
};
//==============================================================================
// Cache contains the cached data for an MGXS object
//==============================================================================
struct CacheDataMG {
int material {-1}; //!< material index
double sqrtkT; //!< last temperature corresponding to t
int t {0}; //!< temperature index
int a {0}; //!< angle index
Direction u; //!< angle that corresponds to a
};
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
@ -231,6 +243,7 @@ private:
vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
vector<ElementMicroXS> photon_xs_; //!< Microscopic photon cross sections
MacroXS macro_xs_; //!< Macroscopic cross sections
CacheDataMG mg_xs_cache_; //!< Multigroup XS cache
int64_t id_; //!< Unique ID
ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.)
@ -349,6 +362,8 @@ public:
ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; }
MacroXS& macro_xs() { return macro_xs_; }
const MacroXS& macro_xs() const { return macro_xs_; }
CacheDataMG& mg_xs_cache() { return mg_xs_cache_; }
const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; }
int64_t& id() { return id_; }
const int64_t& id() const { return id_; }

View file

@ -5,10 +5,6 @@
#include <cstdlib>
#include <sstream>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xadapt.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xsort.hpp"
@ -46,15 +42,6 @@ void Mgxs::init(const std::string& in_name, double in_awr,
n_azi = in_azimuthal.size();
polar = in_polar;
azimuthal = in_azimuthal;
// Set the cross section index cache
#ifdef _OPENMP
int n_threads = omp_get_max_threads();
#else
int n_threads = 1;
#endif
cache.resize(n_threads);
// vector.resize() will value-initialize the members of cache[:]
}
//==============================================================================
@ -425,18 +412,10 @@ void Mgxs::combine(const vector<Mgxs*>& micros, const vector<double>& scalars,
//==============================================================================
double Mgxs::get_xs(
MgxsType xstype, int gin, const int* gout, const double* mu, const int* dg)
double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
const int* dg, int t, int a)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
XsData* xs_t = &xs[cache[tid].t];
int a = cache[tid].a;
#else
XsData* xs_t = &xs[cache[0].t];
int a = cache[0].a;
#endif
XsData* xs_t = &xs[t];
double val;
switch (xstype) {
case MgxsType::TOTAL:
@ -538,19 +517,14 @@ double Mgxs::get_xs(
//==============================================================================
void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed)
void Mgxs::sample_fission_energy(
int gin, int& dg, int& gout, uint64_t* seed, int t, int a)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
XsData* xs_t = &xs[cache[tid].t];
double nu_fission = xs_t->nu_fission(cache[tid].a, gin);
XsData* xs_t = &xs[t];
double nu_fission = xs_t->nu_fission(a, gin);
// Find the probability of having a prompt neutron
double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin);
double prob_prompt = xs_t->prompt_nu_fission(a, gin);
// sample random numbers
double xi_pd = prn(seed) * nu_fission;
@ -566,7 +540,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed)
// sample the outgoing energy group
double prob_gout = 0.;
for (gout = 0; gout < num_groups; ++gout) {
prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout);
prob_gout += xs_t->chi_prompt(a, gin, gout);
if (xi_gout < prob_gout)
break;
}
@ -576,7 +550,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed)
// get the delayed group
for (dg = 0; dg < num_delayed_groups; ++dg) {
prob_prompt += xs_t->delayed_nu_fission(cache[tid].a, dg, gin);
prob_prompt += xs_t->delayed_nu_fission(a, dg, gin);
if (xi_pd < prob_prompt)
break;
}
@ -587,7 +561,7 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed)
// sample the outgoing energy group
double prob_gout = 0.;
for (gout = 0; gout < num_groups; ++gout) {
prob_gout += xs_t->chi_delayed(cache[tid].a, dg, gin, gout);
prob_gout += xs_t->chi_delayed(a, dg, gin, gout);
if (xi_gout < prob_gout)
break;
}
@ -597,35 +571,39 @@ void Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed)
//==============================================================================
void Mgxs::sample_scatter(
int gin, int& gout, double& mu, double& wgt, uint64_t* seed)
int gin, int& gout, double& mu, double& wgt, uint64_t* seed, int t, int a)
{
// This method assumes that the temperature and angle indices are set
// Sample the data
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, seed);
xs[t].scatter[a]->sample(gin, gout, mu, wgt, seed);
}
//==============================================================================
void Mgxs::calculate_xs(Particle& p)
{
// Set our indices
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
set_temperature_index(p.sqrtkT());
set_angle_index(p.u_local());
XsData* xs_t = &xs[cache[tid].t];
p.macro_xs().total = xs_t->total(cache[tid].a, p.g());
p.macro_xs().absorption = xs_t->absorption(cache[tid].a, p.g());
// If the material is different, then we need to do a full lookup
if (p.material() != p.mg_xs_cache().material) {
set_temperature_index(p);
set_angle_index(p);
p.mg_xs_cache().material = p.material();
} else {
// If material is the same, but temperature is different, need to
// find the new temperature index
if (p.sqrtkT() != p.mg_xs_cache().sqrtkT) {
set_temperature_index(p);
}
// If the material is the same, but angle is different, need to
// find the new angle index
if (p.u_local() != p.mg_xs_cache().u) {
set_angle_index(p);
}
}
int temperature = p.mg_xs_cache().t;
int angle = p.mg_xs_cache().a;
p.macro_xs().total = xs[temperature].total(angle, p.g());
p.macro_xs().absorption = xs[temperature].absorption(angle, p.g());
p.macro_xs().nu_fission =
fissionable ? xs_t->nu_fission(cache[tid].a, p.g()) : 0.;
fissionable ? xs[temperature].nu_fission(angle, p.g()) : 0.;
}
//==============================================================================
@ -643,32 +621,26 @@ bool Mgxs::equiv(const Mgxs& that)
//==============================================================================
void Mgxs::set_temperature_index(double sqrtkT)
int Mgxs::get_temperature_index(double sqrtkT) const
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (sqrtkT != cache[tid].sqrtkT) {
cache[tid].t = xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
cache[tid].sqrtkT = sqrtkT;
}
return xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
}
//==============================================================================
void Mgxs::set_angle_index(Direction u)
void Mgxs::set_temperature_index(Particle& p)
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (!is_isotropic && ((u.x != cache[tid].u) || (u.y != cache[tid].v) ||
(u.z != cache[tid].w))) {
p.mg_xs_cache().t = get_temperature_index(p.sqrtkT());
p.mg_xs_cache().sqrtkT = p.sqrtkT();
}
//==============================================================================
int Mgxs::get_angle_index(const Direction& u) const
{
if (is_isotropic) {
return 0;
} else {
// convert direction to polar and azimuthal angles
double my_pol = std::acos(u.z);
double my_azi = std::atan2(u.y, u.x);
@ -679,12 +651,18 @@ void Mgxs::set_angle_index(Direction u)
delta_angle = 2. * PI / n_azi;
int a = std::floor((my_azi + PI) / delta_angle);
cache[tid].a = n_azi * p + a;
return n_azi * p + a;
}
}
// store this direction as the last one used
cache[tid].u = u.x;
cache[tid].v = u.y;
cache[tid].w = u.z;
//==============================================================================
void Mgxs::set_angle_index(Particle& p)
{
// See if we need to find the new index
if (!is_isotropic) {
p.mg_xs_cache().a = get_angle_index(p.u_local());
p.mg_xs_cache().u = p.u_local();
}
}

View file

@ -72,8 +72,8 @@ void sample_reaction(Particle& p)
void scatter(Particle& p)
{
data::mg.macro_xs_[p.material()].sample_scatter(
p.g_last(), p.g(), p.mu(), p.wgt(), p.current_seed());
data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(),
p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);
// Rotate the angle
p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());
@ -145,7 +145,7 @@ void create_fission_sites(Particle& p)
int dg;
int gout;
data::mg.macro_xs_[p.material()].sample_fission_energy(
p.g(), dg, gout, p.current_seed());
p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);
// Store the energy and delayed groups on the fission bank
site.E = gout;

View file

@ -907,7 +907,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
int m;
switch (score_bin) {
// clang-format off
// clang-format off
case N_GAMMA: m = 0; break;
case N_P: m = 1; break;
case N_A: m = 2; break;
@ -1595,10 +1595,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
auto& macro_xs = data::mg.macro_xs_[p.material()];
// Find the temperature and angle indices of interest
macro_xs.set_angle_index(p_u);
int macro_t = p.mg_xs_cache().t;
int macro_a = macro_xs.get_angle_index(p_u);
int nuc_t = 0;
int nuc_a = 0;
if (i_nuclide >= 0) {
nuc_xs.set_temperature_index(p.sqrtkT());
nuc_xs.set_angle_index(p_u);
nuc_t = nuc_xs.get_temperature_index(p.sqrtkT());
nuc_a = nuc_xs.get_angle_index(p_u);
}
for (auto i = 0; i < tally.scores_.size(); ++i) {
@ -1626,12 +1629,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// use the weight of the particle entering the collision as the score
score = flux * p.wgt_last();
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::TOTAL, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score = atom_density * flux * nuc_xs.get_xs(MgxsType::TOTAL, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::TOTAL, p_g, nuc_t, nuc_a);
} else {
score = p.macro_xs().total * flux;
}
@ -1646,17 +1651,21 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// to count 'events' exactly for the inverse velocity
score = flux * p.wgt_last();
if (i_nuclide >= 0) {
score *= nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
score *=
nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a);
} else {
score *= macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
score *=
macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, macro_t, macro_a) /
macro_xs.get_xs(MgxsType::TOTAL, p_g, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score = flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g);
score =
flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g);
score = flux * macro_xs.get_xs(
MgxsType::INVERSE_VELOCITY, p_g, macro_t, macro_a);
}
}
break;
@ -1672,18 +1681,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
if (i_nuclide >= 0) {
score *= atom_density *
nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(),
&p.mu(), nullptr) /
&p.mu(), nullptr, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(),
&p.mu(), nullptr);
&p.mu(), nullptr, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score =
atom_density * flux *
nuc_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, &p.mu(),
nullptr, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(
MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr);
score = flux * macro_xs.get_xs(MgxsType::SCATTER, p_g, nullptr,
&p.mu(), nullptr, macro_t, macro_a);
}
}
break;
@ -1703,16 +1712,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
if (i_nuclide >= 0) {
score *= atom_density *
nuc_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(),
&p.mu(), nullptr) /
&p.mu(), nullptr, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(),
&p.mu(), nullptr);
&p.mu(), nullptr, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score =
atom_density * flux * nuc_xs.get_xs(MgxsType::NU_SCATTER, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::NU_SCATTER, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::NU_SCATTER, p_g);
score =
flux * macro_xs.get_xs(MgxsType::NU_SCATTER, p_g, macro_t, macro_a);
}
}
break;
@ -1732,13 +1742,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = p.wgt_last() * flux;
}
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::ABSORPTION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score =
atom_density * flux * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::ABSORPTION, p_g, nuc_t, nuc_a);
} else {
score = p.macro_xs().absorption * flux;
}
@ -1762,17 +1773,20 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = p.wgt_last() * flux;
}
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
} else {
score *= macro_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *= macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score = atom_density * flux * nuc_xs.get_xs(MgxsType::FISSION, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::FISSION, p_g);
score =
flux * macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
}
}
break;
@ -1793,11 +1807,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// nu-fission
score = wgt_absorb * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *=
atom_density *
nuc_xs.get_xs(MgxsType::NU_FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
} else {
score *= macro_xs.get_xs(MgxsType::NU_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *=
macro_xs.get_xs(MgxsType::NU_FISSION, p_g, macro_t, macro_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
}
} else {
// Skip any non-fission events
@ -1810,16 +1827,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// score.
score = simulation::keff * p.wgt_bank() * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::FISSION, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
}
}
} else {
if (i_nuclide >= 0) {
score =
atom_density * flux * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::NU_FISSION, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::NU_FISSION, p_g);
score =
flux * macro_xs.get_xs(MgxsType::NU_FISSION, p_g, macro_t, macro_a);
}
}
break;
@ -1840,12 +1859,15 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// prompt-nu-fission
score = wgt_absorb * flux;
if (i_nuclide >= 0) {
score *= atom_density *
nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *=
atom_density *
nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
} else {
score *= macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *=
macro_xs.get_xs(
MgxsType::PROMPT_NU_FISSION, p_g, macro_t, macro_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
}
} else {
// Skip any non-fission events
@ -1861,16 +1883,18 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
auto prompt_frac = 1. - n_delayed / static_cast<double>(p.n_bank());
score = simulation::keff * p.wgt_bank() * prompt_frac * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::FISSION, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
}
}
} else {
if (i_nuclide >= 0) {
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g);
nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g);
score = flux * macro_xs.get_xs(
MgxsType::PROMPT_NU_FISSION, p_g, macro_t, macro_a);
}
}
break;
@ -1889,7 +1913,8 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
double abs_xs =
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
if (abs_xs > 0.) {
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -1902,11 +1927,11 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = wgt_absorb * flux;
if (i_nuclide >= 0) {
score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, nuc_t, nuc_a) /
abs_xs;
} else {
score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, macro_t, macro_a) /
abs_xs;
}
score_fission_delayed_dg(
@ -1919,11 +1944,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// delayed-nu-fission xs to the absorption xs
score = wgt_absorb * flux;
if (i_nuclide >= 0) {
score *=
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs;
score *= nuc_xs.get_xs(
MgxsType::DELAYED_NU_FISSION, p_g, nuc_t, nuc_a) /
abs_xs;
} else {
score *=
macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs;
score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
macro_t, macro_a) /
abs_xs;
}
}
}
@ -1948,8 +1975,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = simulation::keff * p.wgt_bank() / p.n_bank() *
p.n_delayed_bank(d - 1) * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::FISSION, p_g);
score *=
atom_density *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
}
score_fission_delayed_dg(
i_tally, d_bin, score, score_index, p.filter_matches());
@ -1962,8 +1991,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score =
simulation::keff * p.wgt_bank() / p.n_bank() * n_delayed * flux;
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::FISSION, p_g);
score *=
atom_density *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
}
}
}
@ -1978,10 +2009,10 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
if (i_nuclide >= 0) {
score = flux * atom_density *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
nullptr, &d);
nullptr, &d, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d);
nullptr, nullptr, &d, macro_t, macro_a);
}
score_fission_delayed_dg(
i_tally, d_bin, score, score_index, p.filter_matches());
@ -1989,10 +2020,12 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
continue;
} else {
if (i_nuclide >= 0) {
score = flux * atom_density *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g);
score =
flux * atom_density *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g);
score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
macro_t, macro_a);
}
}
}
@ -2004,7 +2037,8 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
double abs_xs =
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
if (abs_xs > 0) {
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -2016,16 +2050,16 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
auto d = filt.groups()[d_bin] - 1;
score = wgt_absorb * flux;
if (i_nuclide >= 0) {
score *= nuc_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
score *= nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, nuc_t, nuc_a) *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, nuc_t, nuc_a) /
abs_xs;
} else {
score *= macro_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
score *= macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, macro_t, macro_a) *
macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, macro_t, macro_a) /
abs_xs;
}
score_fission_delayed_dg(
@ -2041,17 +2075,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
if (i_nuclide >= 0) {
score += wgt_absorb * flux *
nuc_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, nuc_t, nuc_a) *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, nuc_t, nuc_a) /
abs_xs;
} else {
score += wgt_absorb * flux *
macro_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, macro_t, macro_a) *
macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d) /
nullptr, nullptr, &d, macro_t, macro_a) /
abs_xs;
}
}
@ -2074,15 +2108,16 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
auto d = bank.delayed_group - 1;
if (d != -1) {
if (i_nuclide >= 0) {
score += simulation::keff * atom_density * bank.wgt * flux *
nuc_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
nuc_xs.get_xs(MgxsType::FISSION, p_g) /
macro_xs.get_xs(MgxsType::FISSION, p_g);
score +=
simulation::keff * atom_density * bank.wgt * flux *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d,
nuc_t, nuc_a) *
nuc_xs.get_xs(MgxsType::FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::FISSION, p_g, macro_t, macro_a);
} else {
score += simulation::keff * bank.wgt * flux *
macro_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d);
macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, macro_t, macro_a);
}
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -2112,17 +2147,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) {
auto d = filt.groups()[d_bin] - 1;
if (i_nuclide >= 0) {
score =
atom_density * flux *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
nuc_xs.get_xs(
MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr,
&d, nuc_t, nuc_a) *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
nullptr, &d, nuc_t, nuc_a);
} else {
score = flux *
macro_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, macro_t, macro_a) *
macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d);
nullptr, nullptr, &d, macro_t, macro_a);
}
score_fission_delayed_dg(
i_tally, d_bin, score, score_index, p.filter_matches());
@ -2132,17 +2167,17 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = 0.;
for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
if (i_nuclide >= 0) {
score +=
atom_density * flux *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
nuc_xs.get_xs(
MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
score += atom_density * flux *
nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, nuc_t, nuc_a) *
nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
nullptr, &d, nuc_t, nuc_a);
} else {
score += flux *
macro_xs.get_xs(
MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) *
macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
nullptr, &d, macro_t, macro_a) *
macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
nullptr, nullptr, &d);
nullptr, nullptr, &d, macro_t, macro_a);
}
}
}
@ -2166,18 +2201,21 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
score = p.wgt_last() * flux;
}
if (i_nuclide >= 0) {
score *= atom_density * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *= atom_density *
nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, nuc_t, nuc_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
} else {
score *= macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
score *=
macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, macro_t, macro_a) /
macro_xs.get_xs(MgxsType::ABSORPTION, p_g, macro_t, macro_a);
}
} else {
if (i_nuclide >= 0) {
score =
atom_density * flux * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g);
score = atom_density * flux *
nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g, nuc_t, nuc_a);
} else {
score = flux * macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g);
score = flux * macro_xs.get_xs(
MgxsType::KAPPA_FISSION, p_g, macro_t, macro_a);
}
}
break;

View file

@ -0,0 +1,53 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0"/>
<macroscopic name="UO2"/>
</material>
<material id="2" name="Water">
<density units="macro" value="1.0"/>
<macroscopic name="LWTR"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" name="fuel inner" region="-1" temperature="600.0" universe="1"/>
<cell id="2" material="1" name="fuel outer" region="1 -2" temperature="294.0" universe="1"/>
<cell id="3" material="2" name="moderator" region="2 3 -4 5 -6" universe="1"/>
<surface coeffs="0.0 0.0 0.25" id="1" name="Fuel IR" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.54" id="2" name="Fuel OR" type="z-cylinder"/>
<surface boundary="reflective" coeffs="-0.63" id="3" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="0.63" id="4" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-0.63" id="5" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="0.63" id="6" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>5</inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="fission">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>
<tallies>
<filter id="1" type="cell">
<bins>1</bins>
</filter>
<filter id="2" type="cell">
<bins>2</bins>
</filter>
<tally id="1" name="inner tally">
<filters>1</filters>
<scores>flux</scores>
</tally>
<tally id="2" name="outer tally">
<filters>2</filters>
<scores>flux</scores>
</tally>
</tallies>
</model>

View file

@ -0,0 +1,8 @@
k-combined:
1.337115E+00 7.221249E-03
tally 1:
2.549066E+01
1.299987E+02
tally 2:
9.276778E+01
1.721791E+03

View file

@ -0,0 +1,166 @@
import os
import numpy as np
import openmc
import openmc.mgxs
from tests.testing_harness import PyAPITestHarness
def create_library():
# Instantiate the energy group data
egroups = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
groups = openmc.mgxs.EnergyGroups(egroups)
# Instantiate the 7-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2', groups, temperatures=[294.0, 600.0])
uo2_xsdata.order = 0
scatter_matrix = np.array([[
[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.4525650, 0.0055664, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]
]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
# Original C5G7 data
uo2_xsdata.set_total([0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, 0.5644058], temperature=294.0)
uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02, 3.0020E-02, 1.1126E-01, 2.8278E-01], temperature=294.0)
uo2_xsdata.set_scatter_matrix(scatter_matrix, temperature=294.0)
uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03, 1.85648E-02, 1.78084E-02, 8.30348E-02, 2.16004E-01], temperature=294.0)
uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02, 4.518301E-02, 4.334208E-02, 2.020901E-01, 5.257105E-01], temperature=294.0)
uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00], temperature=294.0)
# Altered C5G7 data (permuted Chi)
uo2_xsdata.set_total([0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678, 0.5644058], temperature=600.0)
uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02, 3.0020E-02, 1.1126E-01, 2.8278E-01], temperature=600.0)
uo2_xsdata.set_scatter_matrix(scatter_matrix, temperature=600.0)
uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03, 1.85648E-02, 1.78084E-02, 8.30348E-02, 2.16004E-01], temperature=600.0)
uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02, 4.518301E-02, 4.334208E-02, 2.020901E-01, 5.257105E-01], temperature=600.0)
uo2_xsdata.set_chi([4.1176E-01, 5.8791E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00, 0.0000E+00, 0.0000E+00], temperature=600.0)
h2o_xsdata = openmc.XSdata('LWTR', groups)
h2o_xsdata.order = 0
h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
0.718, 1.2544497, 2.650379])
h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02])
scatter_matrix = np.array([[
[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
[0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
[0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
[0.0000000, 0.0000000, 0.0000000, 0.0910284, 0.4155100, 0.0637320, 0.0121390],
[0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]
]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
h2o_xsdata.set_scatter_matrix(scatter_matrix)
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
mg_cross_sections_file.export_to_hdf5()
class MGXSTestHarness(PyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
def test_mg_temperature_multi():
###############################################################################
# Create multigroup data
create_library()
###############################################################################
# Create materials for the problem
# Instantiate some Macroscopic Data
uo2_data = openmc.Macroscopic('UO2')
h2o_data = openmc.Macroscopic('LWTR')
# Instantiate some Materials and register the appropriate Macroscopic objects
uo2 = openmc.Material(name='UO2 fuel')
uo2.set_density('macro', 1.0)
uo2.add_macroscopic(uo2_data)
water = openmc.Material(name='Water')
water.set_density('macro', 1.0)
water.add_macroscopic(h2o_data)
# Instantiate a Materials collection and export to XML
materials = openmc.Materials([uo2, water])
materials.cross_sections = "mgxs.h5"
###############################################################################
# Define problem geometry
# Create a surface for the fuel outer radius
fuel_ir = openmc.ZCylinder(r=0.25, name='Fuel IR')
fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR')
# Create a region represented as the inside of a rectangular prism
pitch = 1.26
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')
# Instantiate Cells
fuel_inner = openmc.Cell(fill=uo2, region=-fuel_ir, name='fuel inner')
fuel_inner.temperature = 600.0
fuel_outer = openmc.Cell(fill=uo2, region=+fuel_ir & -fuel_or, name='fuel outer')
fuel_outer.temperature = 294.0
moderator = openmc.Cell(fill=water, region=+fuel_or & box, name='moderator')
# Create a geometry with the two cells and export to XML
geometry = openmc.Geometry([fuel_inner, fuel_outer, moderator])
###############################################################################
# Define problem settings
# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.energy_mode = "multi-group"
settings.batches = 10
settings.inactive = 5
settings.particles = 1000
# Create an initial uniform spatial source distribution over fissionable zones
lower_left = (-pitch/2, -pitch/2, -1)
upper_right = (pitch/2, pitch/2, 1)
uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True)
settings.source = openmc.IndependentSource(space=uniform_dist)
###############################################################################
# Define tallies
# Instantiate the energy group data
egroups = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6]
inner_filter = openmc.CellFilter(fuel_inner)
outer_filter = openmc.CellFilter(fuel_outer)
energy_filter = openmc.EnergyFilter(egroups)
inner_tally = openmc.Tally(name="inner tally")
inner_tally.filters = [energy_filter]
inner_tally.filters = [inner_filter]
inner_tally.scores = ['flux']
outer_tally = openmc.Tally(name="outer tally")
outer_tally.filters = [energy_filter]
outer_tally.filters = [outer_filter]
outer_tally.scores = ['flux']
# Instantiate a Tallies collection and export to XML
tallies = openmc.Tallies([inner_tally, outer_tally])
# Generate model and run test
model = openmc.Model(geometry, materials, settings, tallies)
harness = MGXSTestHarness('statepoint.10.h5', model=model)
harness.main()