From b5001ee12aafce1f63e4dc8975802922c13892cf Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 21 Jul 2023 02:53:32 -0500 Subject: [PATCH] Multigroup Per-Thread Cache Conversion (#2591) Co-authored-by: Paul Romano --- include/openmc/mgxs.h | 60 ++-- include/openmc/particle_data.h | 15 + src/mgxs.cpp | 142 ++++------ src/physics_mg.cpp | 6 +- src/tallies/tally_scoring.cpp | 262 ++++++++++-------- .../mg_temperature_multi/__init__.py | 0 .../mg_temperature_multi/inputs_true.dat | 53 ++++ .../mg_temperature_multi/results_true.dat | 8 + .../mg_temperature_multi/test.py | 166 +++++++++++ 9 files changed, 487 insertions(+), 225 deletions(-) create mode 100644 tests/regression_tests/mg_temperature_multi/__init__.py create mode 100644 tests/regression_tests/mg_temperature_multi/inputs_true.dat create mode 100644 tests/regression_tests/mg_temperature_multi/results_true.dat create mode 100755 tests/regression_tests/mg_temperature_multi/test.py diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 2da3c83639..d77c34bd5f 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -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 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& get_xsdata() const { return xs; } diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 6d86d528e7..d0813f9c5d 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -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 neutron_xs_; //!< Microscopic neutron cross sections vector 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_; } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 472beda350..9eae5a7da7 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -5,10 +5,6 @@ #include #include -#ifdef _OPENMP -#include -#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& micros, const vector& 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(); } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index ab1fdef028..43e3cfa745 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -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; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 47fe4d1be7..e798161ec2 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -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(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; diff --git a/tests/regression_tests/mg_temperature_multi/__init__.py b/tests/regression_tests/mg_temperature_multi/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/mg_temperature_multi/inputs_true.dat b/tests/regression_tests/mg_temperature_multi/inputs_true.dat new file mode 100644 index 0000000000..b2ad5063c2 --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/inputs_true.dat @@ -0,0 +1,53 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -0.63 -0.63 -1 0.63 0.63 1 + + + multi-group + + + + 1 + + + 2 + + + 1 + flux + + + 2 + flux + + + diff --git a/tests/regression_tests/mg_temperature_multi/results_true.dat b/tests/regression_tests/mg_temperature_multi/results_true.dat new file mode 100644 index 0000000000..80c2f56225 --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/results_true.dat @@ -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 diff --git a/tests/regression_tests/mg_temperature_multi/test.py b/tests/regression_tests/mg_temperature_multi/test.py new file mode 100755 index 0000000000..2090f1b0c2 --- /dev/null +++ b/tests/regression_tests/mg_temperature_multi/test.py @@ -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()