Implementation of migration area score

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

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

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,

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;

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'

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);
 }

src/output.cpp

@@ -609,6 +609,7 @@ const std::unordered_map<int, const char*> 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.

src/reaction.cpp

@@ -187,6 +187,7 @@ std::unordered_map<int, std::string> 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"},

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;

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<int, double>& albedo_map,
   std::unordered_map<int, int>& 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

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,13 +546,18 @@ void Tally::set_scores(const vector<std::string>& 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
     if (!(filt->type() == FilterType::CELL ||
           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<std::string>& 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,13 +1189,22 @@ 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]};
 
     if (tally.active_) {
       model::active_tallies.push_back(i);
       bool mesh_present = (tally.get_filter<MeshFilter>() ||
                            tally.get_filter<MeshMaterialFilter>());
+      if (tally.get_filter<MeshBornFilter>())
+        meshborn_present = true;
+      for (auto score : tally.scores_) {
+        if (score == SCORE_MIGRATION)
+          simulation::migration_present = true;
+      }
       auto time_filter = tally.get_filter<TimeFilter>();
       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()

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

tests/regression_tests/migration/inputs_true.dat

@@ -0,0 +1,29 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="Hydrogen">
+      <density value="1.0" units="g/cm3"/>
+      <nuclide name="H1" ao="1.0"/>
+      <sab name="c_H_in_H2O"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface id="1" type="sphere" boundary="reflective" coeffs="0.0 0.0 0.0 10.0"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>10000</particles>
+    <batches>20</batches>
+    <source type="independent" strength="1.0" particle="neutron"/>
+  </settings>
+  <tallies>
+    <filter id="1" type="energy">
+      <bins>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</bins>
+    </filter>
+    <tally id="1">
+      <filters>1</filters>
+      <scores>migration-area flux total</scores>
+    </tally>
+  </tallies>
+</model>