diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 20a89cea7cb..e69b99920be 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -358,6 +358,16 @@ The following tables show all valid scores: | |kinetics parameters using the iterated fission | | |probability (IFP) method. | +----------------------+---------------------------------------------------+ + |migration-area |The migration-area score is used to evaluate | + | |diffusion coefficients and transport cross sections| + | |in infinite geometry. | + | |The migration area score can only be used with | + | |a ParticleFilter or with an EnergyFilter. | + | |It also cannot be used in the same run as a | + | |MeshBornFilter when non vacuum | + | |boundary conditions are present. | + | |For more information, see Liu_. | + +----------------------+---------------------------------------------------+ .. _usersguide_tally_normalization: @@ -422,3 +432,5 @@ There are several slight variations on this procedure: Note that the only difference between these and the above procedures is in how :math:`H'` is estimated. + +.. _Liu: https://doi.org/10.1016/j.anucene.2017.10.039 diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 0b425a673dc..2d3de989994 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -328,7 +328,8 @@ enum TallyScore { SCORE_PULSE_HEIGHT = -17, // pulse-height SCORE_IFP_TIME_NUM = -18, // IFP lifetime numerator SCORE_IFP_BETA_NUM = -19, // IFP delayed fraction numerator - SCORE_IFP_DENOM = -20 // IFP common denominator + SCORE_IFP_DENOM = -20, // IFP common denominator + SCORE_MIGRATION = -21, // Migration area }; // Global tally parameters diff --git a/include/openmc/position.h b/include/openmc/position.h index 5d291d26b95..20adb8856d3 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -83,6 +83,7 @@ struct Position { return x * other.x + y * other.y + z * other.z; } inline double norm() const { return std::sqrt(x * x + y * y + z * z); } + inline double norm_sq() const { return x * x + y * y + z * z; } inline Position cross(Position other) const { return {y * other.z - z * other.y, z * other.x - x * other.z, diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..5f89fafdf1f 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -32,16 +32,20 @@ extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength extern "C" double - k_abs_tra; //!< sum over batches of k_absorption * k_tracklength + k_abs_tra; //!< sum over batches of k_absorption * k_tracklength +extern bool + migration_present; //! Does an active tally contains a migration score extern double log_spacing; //!< lethargy spacing for energy grid searches extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? -extern "C" int restart_batch; //!< batch at which a restart job resumed -extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? -extern int ssw_current_file; //!< current surface source file -extern "C" int total_gen; //!< total number of generations simulated -extern double total_weight; //!< Total source weight in a batch -extern int64_t work_per_rank; //!< number of particles per MPI rank +extern bool + nonvacuum_boundary_present; //!< Does the geometry contain non-vacuum b.c. +extern "C" int restart_batch; //!< batch at which a restart job resumed +extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? +extern int ssw_current_file; //!< current surface source file +extern "C" int total_gen; //!< total number of generations simulated +extern double total_weight; //!< Total source weight in a batch +extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index c17b16597f9..f049fab524f 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -106,7 +106,7 @@ -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt', -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height', -18: 'ifp-time-numerator', -19: 'ifp-beta-numerator', - -20: 'ifp-denominator', + -20: 'ifp-denominator', -21: 'migration-area', } _ESTIMATORS = { 0: 'analog', 1: 'tracklength', 2: 'collision' diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 5bbda483059..d0b4422a32c 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -7,6 +7,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/random_ray/random_ray.h" +#include "openmc/simulation.h" #include "openmc/surface.h" namespace openmc { @@ -43,6 +44,15 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); + // Handle phantom birth location if migration present + if (simulation::migration_present) { + const auto& r = p.r(); + auto n = surf.normal(r); + n /= n.norm(); + auto& r_born = p.r_born(); + p.r_born() = r_born - 2.0 * (r_born - r).dot(n) * n; + } + p.cross_reflective_bc(surf, u); } @@ -58,6 +68,10 @@ void WhiteBC::handle_particle(Particle& p, const Surface& surf) const // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); + // Handle phantom birth location if migration present + if (simulation::migration_present) + fatal_error("Cannot tally migration area with white boundary conditions."); + p.cross_reflective_bc(surf, u); } @@ -135,6 +149,11 @@ void TranslationalPeriodicBC::handle_particle( // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); + // Handle phantom birth location if migration present + if (simulation::migration_present) { + p.r_born() += translation_; + } + // Pass the new location and surface to the particle. p.cross_periodic_bc(surf, new_r, p.u(), new_surface); } @@ -244,6 +263,18 @@ void RotationalPeriodicBC::handle_particle( // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); + // Handle phantom birth location if migration present + if (simulation::migration_present) { + auto& r_born = p.r_born(); + Position new_r_born; + new_r_born[zero_axis_idx_] = r_born[zero_axis_idx_]; + new_r_born[axis_1_idx_] = + cos_theta * r_born[axis_1_idx_] - sin_theta * r_born[axis_2_idx_]; + new_r_born[axis_2_idx_] = + sin_theta * r_born[axis_1_idx_] + cos_theta * r_born[axis_2_idx_]; + p.r_born() = new_r_born; + } + // Pass the new location, direction, and surface to the particle. p.cross_periodic_bc(surf, new_r, new_u, new_surface); } diff --git a/src/output.cpp b/src/output.cpp index c66c313dce5..592dd23a188 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -609,6 +609,7 @@ const std::unordered_map score_names = { {SCORE_IFP_TIME_NUM, "IFP lifetime numerator"}, {SCORE_IFP_BETA_NUM, "IFP delayed fraction numerator"}, {SCORE_IFP_DENOM, "IFP common denominator"}, + {SCORE_MIGRATION, "Migration Area"}, }; //! Create an ASCII output file showing all tally results. diff --git a/src/reaction.cpp b/src/reaction.cpp index c02e0cc407e..927daa0a4d4 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -187,6 +187,7 @@ std::unordered_map REACTION_NAME_MAP { {SCORE_FLUX, "flux"}, {SCORE_TOTAL, "total"}, {SCORE_SCATTER, "scatter"}, + {SCORE_MIGRATION, "migration-area"}, {SCORE_NU_SCATTER, "nu-scatter"}, {SCORE_ABSORPTION, "absorption"}, {SCORE_FISSION, "fission"}, diff --git a/src/simulation.cpp b/src/simulation.cpp index 4fad196a604..0f5874a69c8 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -309,8 +309,10 @@ double k_col_abs {0.0}; double k_col_tra {0.0}; double k_abs_tra {0.0}; double log_spacing; +bool migration_present {false}; int n_lost_particles {0}; bool need_depletion_rx {false}; +bool nonvacuum_boundary_present {false}; int restart_batch; bool satisfy_triggers {false}; int ssw_current_file; diff --git a/src/surface.cpp b/src/surface.cpp index 81b756deae7..6548d7c5e35 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -17,6 +17,7 @@ #include "openmc/math_functions.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" @@ -1175,6 +1176,8 @@ void read_surfaces(pugi::xml_node node, std::unordered_map& albedo_map, std::unordered_map& periodic_sense_map) { + simulation::nonvacuum_boundary_present = false; + // Count the number of surfaces int n_surfaces = 0; for (pugi::xml_node surf_node : node.children("surface")) { @@ -1245,6 +1248,8 @@ void read_surfaces(pugi::xml_node node, // Check for a periodic surface if (check_for_node(surf_node, "boundary")) { std::string surf_bc = get_node_value(surf_node, "boundary", true, true); + if (surf_bc != "vacuum") + simulation::nonvacuum_boundary_present = true; if (surf_bc == "periodic") { periodic_sense_map[model::surfaces.back()->id_] = 0; // Check for surface albedo. Skip sanity check as it is already done diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 3fe48c1b021..01a1e637591 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -273,6 +273,7 @@ Tally::Tally(pugi::xml_node node) case SCORE_FLUX: case SCORE_TOTAL: case SCORE_SCATTER: + case SCORE_MIGRATION: case SCORE_NU_SCATTER: case SCORE_ABSORPTION: case SCORE_FISSION: @@ -545,6 +546,7 @@ void Tally::set_scores(const vector& scores) bool surface_present = false; bool meshsurface_present = false; bool non_cell_energy_present = false; + bool non_particle_energy_present = false; for (auto i_filt : filters_) { const auto* filt {model::tally_filters[i_filt].get()}; // Checking for only cell and energy filters for pulse-height tally @@ -552,6 +554,10 @@ void Tally::set_scores(const vector& scores) filt->type() == FilterType::ENERGY)) { non_cell_energy_present = true; } + if (!(filt->type() == FilterType::PARTICLE || + filt->type() == FilterType::ENERGY)) { + non_particle_energy_present = true; + } if (filt->type() == FilterType::LEGENDRE) { legendre_present = true; } else if (filt->type() == FilterType::CELLFROM) { @@ -621,6 +627,19 @@ void Tally::set_scores(const vector& scores) estimator_ = TallyEstimator::ANALOG; break; + case SCORE_MIGRATION: + if (estimator_ != TallyEstimator::TRACKLENGTH) + fatal_error( + "Migration-area can only be tallied with a tracklength estimator"); + if (non_particle_energy_present) + fatal_error("Cannot tally migration area with filters other than " + "energy filter and particle filter"); + for (auto i_nuclide : nuclides_) { + if (i_nuclide != NUCLIDE_NONE) + fatal_error("Cannot tally migration area with nuclides."); + } + break; + case SCORE_NU_SCATTER: if (settings::run_CE) { estimator_ = TallyEstimator::ANALOG; @@ -1170,6 +1189,9 @@ void setup_active_tallies() model::active_pulse_height_tallies.clear(); model::time_grid.clear(); + bool meshborn_present = false; + simulation::migration_present = false; + for (auto i = 0; i < model::tallies.size(); ++i) { const auto& tally {*model::tallies[i]}; @@ -1177,6 +1199,12 @@ void setup_active_tallies() model::active_tallies.push_back(i); bool mesh_present = (tally.get_filter() || tally.get_filter()); + if (tally.get_filter()) + meshborn_present = true; + for (auto score : tally.scores_) { + if (score == SCORE_MIGRATION) + simulation::migration_present = true; + } auto time_filter = tally.get_filter(); switch (tally.type_) { @@ -1212,6 +1240,10 @@ void setup_active_tallies() } } } + if (meshborn_present && simulation::migration_present && + simulation::nonvacuum_boundary_present) + fatal_error("Cannot score migration-area in the same simulation as a " + "MeshBorn filter and a non vacuum b.c."); } void free_memory_tally() diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 4210b034a96..af2ae81d290 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -670,6 +670,12 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, } break; + case SCORE_MIGRATION: + score = 1.0 / 6.0 * p.wgt() * + (p.r() - p.r_last()) + .dot((p.r() - p.r_born()) + (p.r_last() - p.r_born())); + break; + case SCORE_NU_FISSION: if (p.macro_xs().fission == 0) continue; @@ -1767,6 +1773,12 @@ void score_general_mg(Particle& p, int i_tally, int start_index, } break; + case SCORE_MIGRATION: + score = 1.0 / 6.0 * p.wgt() * + (p.r() - p.r_last()) + .dot((p.r() - p.r_born()) + (p.r_last() - p.r_born())); + break; + case SCORE_NU_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter diff --git a/tests/regression_tests/migration/__init__.py b/tests/regression_tests/migration/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/migration/inputs_true.dat b/tests/regression_tests/migration/inputs_true.dat new file mode 100644 index 00000000000..2ecc3d2a02f --- /dev/null +++ b/tests/regression_tests/migration/inputs_true.dat @@ -0,0 +1,29 @@ + + + + + + + + + + + + + + + fixed source + 10000 + 20 + + + + + 0.0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.042 0.05 0.058 0.067 0.08 0.1 0.14 0.18 0.22 0.25 0.28 0.3 0.32 0.35 0.4 0.5 0.625 0.78 0.85 0.91 0.95 0.972 0.996 1.02 1.045 1.071 1.097 1.123 1.15 1.3 1.5 1.855 2.1 2.6 3.3 4.0 9.877 15.968 27.7 48.052 75.501 148.73 367.26 906.9 1425.1 2239.5 3519.1 5530.0 9118.0 15030.0 24780.0 40850.0 67340.0 111000.0 183000.0 302500.0 500000.0 821000.0 1353000.0 2231000.0 3679000.0 6065500.0 20000000.0 + + + 1 + migration-area flux total + + + diff --git a/tests/regression_tests/migration/results_true.dat b/tests/regression_tests/migration/results_true.dat new file mode 100644 index 00000000000..825d7893b06 --- /dev/null +++ b/tests/regression_tests/migration/results_true.dat @@ -0,0 +1,421 @@ +tally 1: +1.151715E-02 +7.764946E-06 +1.885172E+00 +1.777098E-01 +1.033112E+02 +5.337084E+02 +3.983862E-02 +8.306135E-05 +4.812785E+00 +1.158257E+00 +2.164162E+02 +2.342026E+03 +6.102772E-02 +1.933110E-04 +6.390797E+00 +2.042296E+00 +2.539320E+02 +3.224358E+03 +7.454985E-02 +2.832220E-04 +7.506089E+00 +2.817389E+00 +2.692124E+02 +3.624175E+03 +9.519027E-02 +4.579944E-04 +7.933918E+00 +3.147711E+00 +2.603718E+02 +3.390061E+03 +1.084204E-01 +5.959081E-04 +8.066798E+00 +3.253956E+00 +2.451357E+02 +3.004846E+03 +1.170748E-01 +6.927723E-04 +7.857292E+00 +3.087114E+00 +2.232923E+02 +2.493184E+03 +1.668098E-01 +1.401994E-03 +1.030060E+01 +5.305467E+00 +2.738330E+02 +3.749472E+03 +1.914057E-01 +1.841487E-03 +1.053553E+01 +5.550439E+00 +2.619708E+02 +3.431791E+03 +1.747476E-01 +1.536344E-03 +9.019503E+00 +4.067920E+00 +2.125963E+02 +2.260052E+03 +1.706312E-01 +1.467254E-03 +8.477252E+00 +3.593542E+00 +1.914483E+02 +1.832801E+03 +2.088285E-01 +2.199507E-03 +9.252138E+00 +4.280461E+00 +2.004146E+02 +2.008469E+03 +2.375246E-01 +2.835795E-03 +9.804908E+00 +4.807318E+00 +2.017756E+02 +2.035882E+03 +2.399958E-01 +2.893945E-03 +8.651606E+00 +3.742882E+00 +1.660699E+02 +1.379097E+03 +8.764552E-02 +3.897486E-04 +2.771932E+00 +3.842349E-01 +4.931669E+01 +1.216241E+02 +4.015863E-02 +8.219874E-05 +1.043561E+00 +5.447041E-02 +1.752695E+01 +1.536519E+01 +1.629873E-02 +1.432929E-05 +3.954570E-01 +7.824049E-03 +6.390631E+00 +2.043242E+00 +1.392893E-02 +1.021488E-05 +2.770772E-01 +3.845190E-03 +4.360695E+00 +9.524174E-01 +8.125582E-03 +3.478406E-06 +1.537727E-01 +1.185045E-03 +2.372366E+00 +2.820595E-01 +8.004140E-03 +3.297562E-06 +1.404136E-01 +9.888175E-04 +2.134762E+00 +2.285577E-01 +1.048598E-02 +5.754071E-06 +1.846759E-01 +1.706879E-03 +2.760633E+00 +3.814141E-01 +1.598730E-02 +1.339635E-05 +2.710495E-01 +3.679302E-03 +3.951106E+00 +7.818129E-01 +2.983250E-02 +4.514200E-05 +4.452826E-01 +9.923714E-03 +6.243390E+00 +1.950934E+00 +3.087791E-02 +4.885840E-05 +4.330658E-01 +9.382965E-03 +5.854196E+00 +1.714615E+00 +3.197897E-02 +5.185329E-05 +4.196693E-01 +8.815073E-03 +5.552245E+00 +1.542943E+00 +1.222073E-02 +7.772441E-06 +1.602311E-01 +1.286219E-03 +2.097335E+00 +2.203719E-01 +1.049044E-02 +5.848900E-06 +1.270539E-01 +8.085071E-04 +1.655124E+00 +1.372048E-01 +5.409500E-03 +1.611747E-06 +7.662452E-02 +2.949919E-04 +9.947459E-01 +4.971637E-02 +2.954315E-03 +4.979504E-07 +4.135712E-02 +8.594375E-05 +5.358078E-01 +1.442556E-02 +3.475106E-03 +6.835279E-07 +4.413458E-02 +9.833210E-05 +5.709820E-01 +1.645823E-02 +3.687034E-03 +7.651751E-07 +4.276488E-02 +9.203002E-05 +5.524587E-01 +1.535869E-02 +2.973344E-03 +5.655796E-07 +4.305110E-02 +9.362585E-05 +5.553167E-01 +1.557791E-02 +3.438046E-03 +6.888869E-07 +4.453476E-02 +9.999720E-05 +5.736155E-01 +1.658943E-02 +3.534723E-03 +7.437425E-07 +4.353717E-02 +9.564247E-05 +5.600259E-01 +1.582513E-02 +3.336017E-03 +6.338724E-07 +4.281204E-02 +9.252136E-05 +5.499627E-01 +1.526782E-02 +3.132302E-03 +5.615805E-07 +4.117330E-02 +8.567774E-05 +5.282056E-01 +1.410080E-02 +1.611275E-02 +1.359495E-05 +2.187754E-01 +2.396084E-03 +2.796225E+00 +3.914248E-01 +2.084727E-02 +2.198299E-05 +2.542216E-01 +3.234727E-03 +3.230357E+00 +5.222926E-01 +2.832940E-02 +4.066302E-05 +3.631286E-01 +6.603904E-03 +4.584601E+00 +1.052648E+00 +1.742385E-02 +1.562468E-05 +2.121201E-01 +2.254130E-03 +2.664219E+00 +3.555945E-01 +2.823366E-02 +4.044399E-05 +3.568934E-01 +6.371694E-03 +4.463419E+00 +9.965817E-01 +3.306521E-02 +5.574268E-05 +3.960799E-01 +7.849005E-03 +4.930000E+00 +1.216025E+00 +2.583372E-02 +3.425765E-05 +3.214348E-01 +5.169414E-03 +3.986209E+00 +7.950156E-01 +1.288130E-01 +8.347117E-04 +1.487536E+00 +1.106615E-01 +1.824114E+01 +1.664045E+01 +6.724831E-02 +2.269465E-04 +7.856337E-01 +3.087087E-02 +9.609752E+00 +4.618841E+00 +7.631762E-02 +2.922324E-04 +9.058144E-01 +4.103555E-02 +1.107279E+01 +6.131912E+00 +7.669017E-02 +2.954654E-04 +9.047361E-01 +4.094391E-02 +1.105412E+01 +6.112144E+00 +6.188300E-02 +1.937968E-04 +7.422126E-01 +2.755300E-02 +9.064836E+00 +4.109910E+00 +9.483947E-02 +4.524101E-04 +1.110980E+00 +6.173095E-02 +1.356214E+01 +9.199129E+00 +1.313218E-01 +8.666723E-04 +1.486576E+00 +1.105135E-01 +1.812858E+01 +1.643496E+01 +1.370542E-01 +9.442553E-04 +1.496134E+00 +1.119550E-01 +1.820094E+01 +1.656877E+01 +7.434513E-02 +2.778881E-04 +7.458280E-01 +2.782318E-02 +9.040309E+00 +4.087858E+00 +7.539331E-02 +2.850294E-04 +7.457057E-01 +2.781265E-02 +9.000199E+00 +4.051460E+00 +7.784473E-02 +3.051218E-04 +7.596350E-01 +2.887538E-02 +9.107637E+00 +4.150777E+00 +8.326257E-02 +3.473454E-04 +7.528319E-01 +2.835072E-02 +8.933289E+00 +3.992002E+00 +9.989220E-02 +5.007887E-04 +8.565505E-01 +3.670562E-02 +9.991295E+00 +4.994253E+00 +1.097272E-01 +6.034104E-04 +8.781243E-01 +3.856616E-02 +9.957477E+00 +4.958973E+00 +1.267186E-01 +8.044676E-04 +9.225617E-01 +4.257058E-02 +1.001088E+01 +5.012587E+00 +1.484297E-01 +1.104237E-03 +9.826252E-01 +4.829042E-02 +9.975955E+00 +4.977335E+00 +1.876936E-01 +1.764468E-03 +1.092320E+00 +5.968028E-02 +1.002380E+01 +5.025684E+00 +2.428358E-01 +2.950948E-03 +1.257521E+00 +7.908212E-02 +1.003826E+01 +5.039217E+00 +3.227289E-01 +5.213680E-03 +1.517685E+00 +1.152121E-01 +1.012758E+01 +5.130297E+00 +4.609490E-01 +1.063557E-02 +1.924087E+00 +1.852168E-01 +1.031358E+01 +5.321697E+00 +6.510120E-01 +2.121435E-02 +2.476179E+00 +3.067256E-01 +1.033012E+01 +5.338175E+00 +9.296729E-01 +4.325819E-02 +3.248571E+00 +5.278302E-01 +1.042151E+01 +5.432160E+00 +1.279432E+00 +8.193963E-02 +4.207982E+00 +8.856296E-01 +1.037174E+01 +5.380346E+00 +1.592958E+00 +1.270860E-01 +4.899777E+00 +1.201256E+00 +9.212557E+00 +4.246567E+00 +1.695755E+00 +1.440814E-01 +4.697417E+00 +1.104031E+00 +6.645772E+00 +2.209786E+00 +1.246831E+00 +7.812329E-02 +2.927301E+00 +4.290505E-01 +3.035756E+00 +4.613988E-01 +4.570605E-01 +1.061146E-02 +8.400634E-01 +3.542841E-02 +6.054760E-01 +1.840507E-02 diff --git a/tests/regression_tests/migration/test.py b/tests/regression_tests/migration/test.py new file mode 100644 index 00000000000..f84af544a2a --- /dev/null +++ b/tests/regression_tests/migration/test.py @@ -0,0 +1,44 @@ +"""Test migration area calculation.""" + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + +@pytest.fixture() +def model(): + model = openmc.Model() + + # Material + material = openmc.Material(name="Hydrogen") + material.add_nuclide("H1", 1.0) + material.set_density('g/cm3', 1.0) + material.add_s_alpha_beta('c_H_in_H2O') + + # Geometry + radius = 10.0 + sphere = openmc.Sphere(r=radius, boundary_type="reflective") + cell = openmc.Cell(region=-sphere, fill=material) + model.geometry = openmc.Geometry([cell]) + + # Settings + model.settings.particles = 10000 + model.settings.batches = 20 + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.IndependentSource() + + + # Tally + groups = openmc.mgxs.EnergyGroups('CASMO-70') + energy_filter = openmc.EnergyFilter(groups.group_edges) + tally = openmc.Tally() + tally.scores = ["migration-area", "flux", "total"] + tally.filters = [energy_filter] + model.tallies = [tally] + + return model + + +def test_migration_area(model): + harness = PyAPITestHarness("statepoint.20.h5", model=model) + harness.main() diff --git a/tests/unit_tests/test_migration.py b/tests/unit_tests/test_migration.py new file mode 100644 index 00000000000..a0aa2203616 --- /dev/null +++ b/tests/unit_tests/test_migration.py @@ -0,0 +1,53 @@ +import openmc +import pytest +import numpy as np + +from tests import cdtemp + +def sphere_model(radius, boundary_type): + + model = openmc.Model() + + # Material + material = openmc.Material(name="Hydrogen") + material.add_nuclide("H1", 1.0) + material.set_density('g/cm3', 1.0) + material.add_s_alpha_beta('c_H_in_H2O') + + # Geometry + sphere = openmc.Sphere(r=radius, boundary_type=boundary_type) + cell = openmc.Cell(region=-sphere, fill=material) + model.geometry = openmc.Geometry([cell]) + + # Settings + model.settings.particles = 1000 + model.settings.batches = 5 + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.IndependentSource() + + + # Tally + tally = openmc.Tally() + tally.scores = ["migration-area"] + model.tallies = [tally] + + return model, tally + +def test_reflection_is_equivalent_to_large_model(run_in_tmpdir): + openmc.reset_auto_ids() + refl_model, tally1 = sphere_model(2.5, "reflective") + large_model, tally2 = sphere_model(100, "vacuum") + + with cdtemp(): + refl_model.run(apply_tally_results=True) + mean1, std1 = (tally1.mean.squeeze(), + tally1.std_dev.squeeze()) + + with cdtemp(): + large_model.run(apply_tally_results=True) + mean2, std2 = (tally2.mean.squeeze(), + tally2.std_dev.squeeze()) + std = (std1**2+std2**2)**0.5 + diff = np.abs(mean2-mean1) + + assert diff <= 3*std