Moved flux derivatives to the particle. Diff tallies are still not getting right answer, but they look closer...

include/openmc/particle.h

@@ -313,6 +313,8 @@ public:
   std::vector<Particle::Bank> secondary_bank_;
 
   int64_t current_work_; // current work index
+
+  std::vector<double> flux_derivs_;  // Derivatives of the current particle's weight
 };
 
 } // namespace openmc

include/openmc/tallies/derivative.h

@@ -19,7 +19,6 @@ struct TallyDerivative {
   int variable;  //!< Independent variable (like temperature)
   int diff_material;  //!< Material this derivative is applied to
   int diff_nuclide;  //!< Nuclide this material is applied to
-  double flux_deriv;  //!< Derivative of the current particle's weight
 
   TallyDerivative() {}
   explicit TallyDerivative(pugi::xml_node node);
@@ -46,16 +45,16 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
 //! further tallies are scored.
 //
 //! \param p The particle being tracked
-void score_collision_derivative(const Particle* p);
+void score_collision_derivative(Particle* p);
 
 //! Adjust diff tally flux derivatives for a particle tracking event.
 //
 //! \param p The particle being tracked
 //! \param distance The distance in [cm] traveled by the particle
-void score_track_derivative(const Particle* p, double distance);
+void score_track_derivative(Particle* p, double distance);
 
 //! Set the flux derivatives on differential tallies to zero.
-void zero_flux_derivs();
+void zero_flux_derivs(std::vector<double> v);
 
 } // namespace openmc
 
@@ -63,15 +62,10 @@ void zero_flux_derivs();
 // Global variables
 //==============================================================================
 
-// Explicit vector template specialization declaration of threadprivate variable
-// outside of the openmc namespace for the picky Intel compiler.
-//extern template class std::vector<openmc::TallyDerivative>;
-
 namespace openmc {
 
 namespace model {
 extern std::vector<TallyDerivative> tally_derivs;
-//#pragma omp threadprivate(tally_derivs)
 extern std::unordered_map<int, int> tally_deriv_map;
 } // namespace model
 

src/particle.cpp

@@ -109,6 +109,7 @@ Particle::from_source(const Bank* src)
   material_ = C_NONE;
   n_collision_ = 0;
   fission_ = false;
+  std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0);
 
   // Copy attributes from source bank site
   type_ = src->particle;
@@ -154,7 +155,8 @@ Particle::transport()
   if (write_track_) add_particle_track();
 
   // Every particle starts with no accumulated flux derivative.
-  if (!model::active_tallies.empty()) zero_flux_derivs();
+  if (!model::active_tallies.empty())
+    flux_derivs_.resize(model::tally_derivs.size(), 0.0);
 
   while (true) {
     // Set the random number stream

src/simulation.cpp

@@ -1280,35 +1280,29 @@ void initialize_history(Particle* p, int64_t index_source)
     }
   }
 
-// Display message if high verbosity or trace is on
-   if (settings::verbosity >= 9 || simulation::trace) {
-      write_message("Simulating Particle " + std::to_string(p->id_));
-   }
+  // Display message if high verbosity or trace is on
+  if (settings::verbosity >= 9 || simulation::trace) {
+    write_message("Simulating Particle " + std::to_string(p->id_));
+  }
 
-    // // Initialize number of events to zero
-   // int n_event = 0;
+  // Add paricle's starting weight to count for normalizing tallies later
+  #pragma omp atomic
+  simulation::total_weight += p->wgt_;
 
-    // Add paricle's starting weight to count for normalizing tallies later
-   #pragma omp atomic
-   simulation::total_weight += p->wgt_;
+  // Force calculation of cross-sections by setting last energy to zero
+  if (settings::run_CE) {
+    for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
+  }
 
-    // Force calculation of cross-sections by setting last energy to zero
-   if (settings::run_CE) {
-     for (auto& micro : p->neutron_xs_) micro.last_E = 0.0;
-   }
+  // Prepare to write out particle track.
+  if (p->write_track_) add_particle_track();
 
-    // Prepare to write out particle track.
-   if (p->write_track_) add_particle_track();
-
-    // Every particle starts with no accumulated flux derivative.
-   if (!model::active_tallies.empty()) zero_flux_derivs();
-
-   /*
-   std::cout << "Initialized particle " << particle_seed << " with E = " << p->E_ << " and Position {" <<
-	   p->r().x << ", " <<
-	   p->r().y << ", " <<
-	   p->r().z << "}" << std::endl;
-	   */
+  // Every particle starts with no accumulated flux derivative.
+  if (!model::active_tallies.empty())
+  {
+    p->flux_derivs_.resize(model::tally_derivs.size(), 0.0);
+    std::fill(p->flux_derivs_.begin(), p->flux_derivs_.end(), 0.0);
+  }
 }
 
 int overall_generation()

src/tallies/derivative.cpp

@@ -81,13 +81,9 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
 void
 read_tally_derivatives(pugi::xml_node node)
 {
-  // Populate the derivatives array.  This must be done in parallel because
-  // the derivatives are threadprivate.
-  //#pragma omp parallel
-  {
-    for (auto deriv_node : node.children("derivative"))
-      model::tally_derivs.emplace_back(deriv_node);
-  }
+  // Populate the derivatives array.
+  for (auto deriv_node : node.children("derivative"))
+    model::tally_derivs.emplace_back(deriv_node);
 
   // Fill the derivative map.
   for (auto i = 0; i < model::tally_derivs.size(); ++i) {
@@ -119,8 +115,8 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
   // where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the
   // perturbated variable.
 
-  const auto& deriv {model::tally_derivs[tally.deriv_]};
-  auto flux_deriv = deriv.flux_deriv;
+  const auto&   deriv      {model::tally_derivs[tally.deriv_]};
+  const double& flux_deriv {p->flux_derivs_[tally.deriv_]};
 
   // Handle special cases where we know that d_c/d_p must be zero.
   if (score_bin == SCORE_FLUX) {
@@ -567,13 +563,15 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
 }
 
 void
-score_track_derivative(const Particle* p, double distance)
+score_track_derivative(Particle* p, double distance)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material_ == MATERIAL_VOID) return;
   const Material& material {*model::materials[p->material_]};
-
-  for (auto& deriv : model::tally_derivs) {
+  
+  for (int idx = 0; idx < model::tally_derivs.size(); idx++) {
+    const auto& deriv {model::tally_derivs[idx]};
+    double& flux_deriv = p->flux_derivs_[idx];
     if (deriv.diff_material != material.id_) continue;
 
     switch (deriv.variable) {
@@ -582,7 +580,7 @@ score_track_derivative(const Particle* p, double distance)
       // phi is proportional to e^(-Sigma_tot * dist)
       // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
       // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
-      deriv.flux_deriv -= distance * p->macro_xs_.total
+      flux_deriv -= distance * p->macro_xs_.total
         / material.density_gpcc_;
       break;
 
@@ -590,7 +588,7 @@ score_track_derivative(const Particle* p, double distance)
       // phi is proportional to e^(-Sigma_tot * dist)
       // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
       // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
-      deriv.flux_deriv -= distance
+      flux_deriv -= distance
         * p->neutron_xs_[deriv.diff_nuclide].total;
       break;
 
@@ -604,7 +602,7 @@ score_track_derivative(const Particle* p, double distance)
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_);
-          deriv.flux_deriv -= distance * (dsig_s + dsig_a)
+          flux_deriv -= distance * (dsig_s + dsig_a)
             * material.atom_density_(i);
         }
       }
@@ -613,14 +611,18 @@ score_track_derivative(const Particle* p, double distance)
   }
 }
 
-void score_collision_derivative(const Particle* p)
+void score_collision_derivative(Particle* p)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material_ == MATERIAL_VOID) return;
 
   const Material& material {*model::materials[p->material_]};
 
-  for (auto& deriv : model::tally_derivs) {
+  //for (auto& deriv : model::tally_derivs) {
+  for (int idx = 0; idx < model::tally_derivs.size(); idx++) {
+    const auto&  deriv = model::tally_derivs[idx];
+    double& flux_deriv = p->flux_derivs_[idx];
+
     if (deriv.diff_material != material.id_) continue;
 
     switch (deriv.variable) {
@@ -629,7 +631,7 @@ void score_collision_derivative(const Particle* p)
       // phi is proportional to Sigma_s
       // (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
       // (1 / phi) * (d_phi / d_rho) = 1 / rho
-      deriv.flux_deriv += 1. / material.density_gpcc_;
+      flux_deriv += 1. / material.density_gpcc_;
       break;
 
     case DIFF_NUCLIDE_DENSITY:
@@ -650,7 +652,7 @@ void score_collision_derivative(const Particle* p)
       // (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s
       // (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s
       // (1 / phi) * (d_phi / d_N) = 1 / N
-      deriv.flux_deriv += 1. / material.atom_density_(i);
+      flux_deriv += 1. / material.atom_density_(i);
       break;
 
     case DIFF_TEMPERATURE:
@@ -665,7 +667,7 @@ void score_collision_derivative(const Particle* p)
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
-          deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
+          flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
           // Note that this is an approximation!  The real scattering cross
           // section is
           // Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u).
@@ -680,9 +682,9 @@ void score_collision_derivative(const Particle* p)
   }
 }
 
-void zero_flux_derivs()
+void zero_flux_derivs(std::vector<double> v)
 {
-  for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.;
+  std::fill(v.begin(), v.end(), 0);
 }
 
 }// namespace openmc

src/tallies/tally.cpp

@@ -1042,10 +1042,7 @@ setup_active_tallies()
 void
 free_memory_tally()
 {
-  //#pragma omp parallel
-  {
-    model::tally_derivs.clear();
-  }
+  model::tally_derivs.clear();
   model::tally_deriv_map.clear();
 
   model::tally_filters.clear();