From fe77e3fdbb6d0b582915207b7f2213b598e36d38 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 16 Feb 2026 20:30:04 +0200 Subject: [PATCH 01/10] initial implementation of migration area score --- include/openmc/constants.h | 3 +- include/openmc/position.h | 1 + include/openmc/simulation.h | 18 +- openmc/lib/tally.py | 2 +- src/boundary_condition.cpp | 31 ++ src/output.cpp | 1 + src/reaction.cpp | 1 + src/simulation.cpp | 2 + src/surface.cpp | 5 + src/tallies/tally.cpp | 32 ++ src/tallies/tally_scoring.cpp | 12 + tests/regression_tests/migration/__init__.py | 0 .../migration/inputs_true.dat | 29 ++ .../migration/results_true.dat | 421 ++++++++++++++++++ tests/regression_tests/migration/test.py | 44 ++ 15 files changed, 593 insertions(+), 9 deletions(-) create mode 100644 tests/regression_tests/migration/__init__.py create mode 100644 tests/regression_tests/migration/inputs_true.dat create mode 100644 tests/regression_tests/migration/results_true.dat create mode 100644 tests/regression_tests/migration/test.py diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 38238598e..bd13aa794 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -327,7 +327,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 5d291d26b..20adb8856 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 9a6cf1b21..5f89fafdf 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 c17b16597..f049fab52 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -106,7 +106,7 @@ _SCORES = { -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 5bbda4830..7bdf072ea 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) { + 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 ae2daaffc..3d448803c 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -598,6 +598,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 9ac5a1f52..5043c1404 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 18e40a8bc..627a5cbe7 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -308,8 +308,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 81b756dea..6548d7c5e 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 8cdf3a905..649818386 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -270,6 +270,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: @@ -540,6 +541,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 @@ -547,6 +549,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) { @@ -601,6 +607,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 tallies with 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; @@ -1150,6 +1169,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]}; @@ -1157,6 +1179,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_) { @@ -1192,6 +1220,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 ccc70fc8f..a4cf0f497 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -668,6 +668,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_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq()); + break; + case SCORE_NU_FISSION: if (p.macro_xs().fission == 0) continue; @@ -1765,6 +1771,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_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq()); + 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 000000000..e69de29bb diff --git a/tests/regression_tests/migration/inputs_true.dat b/tests/regression_tests/migration/inputs_true.dat new file mode 100644 index 000000000..2ecc3d2a0 --- /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 000000000..825d7893b --- /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 000000000..f84af544a --- /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() From 566c383f8650f9dcf530763552e8a8fb8094cc9a Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 16 Feb 2026 21:10:58 +0200 Subject: [PATCH 02/10] add score to docs --- docs/source/usersguide/tallies.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index e3b4e508b..b07f5a4ad 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -340,6 +340,15 @@ 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 for non vacuum boundary conditions. | + | |For more information, see Liu_. | + +----------------------+---------------------------------------------------+ .. _usersguide_tally_normalization: @@ -404,3 +413,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 From c81063b52019034eccee4dc3cba33d7ea66b7671 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 10 Mar 2026 21:28:05 +0200 Subject: [PATCH 03/10] Update src/tallies/tally.cpp Co-authored-by: Guillaume Giudicelli --- src/tallies/tally.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 649818386..3a24664f9 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -610,7 +610,7 @@ void Tally::set_scores(const vector& scores) case SCORE_MIGRATION: if (estimator_ != TallyEstimator::TRACKLENGTH) fatal_error( - "Migration-area can only be tallies with tracklength estimator"); + "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"); From d749e0cbade2bd0ed658bc668b66f2beee8ec4c8 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 10 Mar 2026 21:28:44 +0200 Subject: [PATCH 04/10] Update docs/source/usersguide/tallies.rst Co-authored-by: Guillaume Giudicelli --- docs/source/usersguide/tallies.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index b07f5a4ad..35bc231ec 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -346,7 +346,8 @@ The following tables show all valid scores: | |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 for non vacuum boundary conditions. | + | |MeshBornFilter when non vacuum | + | |boundary conditions are present. | | |For more information, see Liu_. | +----------------------+---------------------------------------------------+ From 2ca33428c2844b22b111c1e5b628b04f6594dc86 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 10 Mar 2026 21:51:13 +0200 Subject: [PATCH 05/10] Update src/boundary_condition.cpp Co-authored-by: Guillaume Giudicelli --- src/boundary_condition.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 7bdf072ea..99448a8a3 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -46,10 +46,10 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const // Handle phantom birth location if migration present if (simulation::migration_present) { - auto r = p.r(); + const auto & r = p.r(); auto n = surf.normal(r); n /= n.norm(); - auto r_born = p.r_born(); + auto & r_born = p.r_born(); p.r_born() = r_born - 2.0 * (r_born - r).dot(n) * n; } From d35a8f90aacdf97306d14d709948ac0008fa922d Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Tue, 10 Mar 2026 21:51:33 +0200 Subject: [PATCH 06/10] Update src/boundary_condition.cpp Co-authored-by: Guillaume Giudicelli --- src/boundary_condition.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 99448a8a3..202b3b5b9 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -265,7 +265,7 @@ void RotationalPeriodicBC::handle_particle( // Handle phantom birth location if migration present if (simulation::migration_present) { - auto r_born = p.r_born(); + 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_] = From e43e6e7bcafe6e3023f8e0c5427a61d46a45ec0e Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 10 Mar 2026 21:52:46 +0200 Subject: [PATCH 07/10] apply clang-format --- src/boundary_condition.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 202b3b5b9..d0b4422a3 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -46,10 +46,10 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const // Handle phantom birth location if migration present if (simulation::migration_present) { - const auto & r = p.r(); + const auto& r = p.r(); auto n = surf.normal(r); n /= n.norm(); - auto & r_born = p.r_born(); + auto& r_born = p.r_born(); p.r_born() = r_born - 2.0 * (r_born - r).dot(n) * n; } @@ -265,7 +265,7 @@ void RotationalPeriodicBC::handle_particle( // Handle phantom birth location if migration present if (simulation::migration_present) { - auto & r_born = p.r_born(); + 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_] = From 86a7b8f0cedee531b9f834efc0c34cdb7cb35b80 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 10 Mar 2026 23:34:20 +0200 Subject: [PATCH 08/10] added test and improved numerics --- src/tallies/tally_scoring.cpp | 16 ++++----- tests/unit_tests/test_migration.py | 53 ++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 8 deletions(-) create mode 100644 tests/unit_tests/test_migration.py diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index a4cf0f497..d7c2018a9 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -669,9 +669,9 @@ 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_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq()); + 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: @@ -1772,9 +1772,9 @@ 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_born()).norm_sq() - (p.r_last() - p.r_born()).norm_sq()); + 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: @@ -2646,7 +2646,7 @@ void score_surface_tally(Particle& p, const vector& tallies) // for a further scoring function. double score = current * filter_weight; for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { #pragma omp atomic tally.results_(filter_index, score_index, TallyResult::VALUE) += score; } @@ -2721,7 +2721,7 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) // Loop over scores. for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { #pragma omp atomic tally.results_(filter_index, score_index, TallyResult::VALUE) += filter_weight; diff --git a/tests/unit_tests/test_migration.py b/tests/unit_tests/test_migration.py new file mode 100644 index 000000000..a0aa22036 --- /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 From 148c914c55c792318537077eba6b416d3242da8d Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 10 Mar 2026 23:56:26 +0200 Subject: [PATCH 09/10] ran clang format --- src/tallies/tally_scoring.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index dd7bd1c8d..b010e202e 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2702,7 +2702,7 @@ void score_surface_tally( // Loop over scores. for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { auto score_bin = tally.scores_[score_index]; double score; if (score_bin == SCORE_CURRENT) { @@ -2781,7 +2781,7 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) // Loop over scores. for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { #pragma omp atomic tally.results_(filter_index, score_index, TallyResult::VALUE) += filter_weight; From aa6adfa7b8528beb1cdc8cc0913d5e5a39ef02ad Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Mar 2026 00:00:58 +0200 Subject: [PATCH 10/10] ran clang format --- src/tallies/tally_scoring.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index b010e202e..af2ae81d2 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2648,7 +2648,7 @@ void score_meshsurface_tally(Particle& p, const vector& tallies) // for a further scoring function. double score = current * filter_weight; for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { #pragma omp atomic tally.results_(filter_index, score_index, TallyResult::VALUE) += score; } @@ -2702,7 +2702,7 @@ void score_surface_tally( // Loop over scores. for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { auto score_bin = tally.scores_[score_index]; double score; if (score_bin == SCORE_CURRENT) { @@ -2781,7 +2781,7 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) // Loop over scores. for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + ++score_index) { #pragma omp atomic tally.results_(filter_index, score_index, TallyResult::VALUE) += filter_weight;