SA before merge upstream

Author: Ebiwonjumi

include/openmc/particle_data.h

@@ -48,7 +48,6 @@ struct SourceSite {
   double time {0.0};
   double wgt {1.0};
   int delayed_group {0};
-  int surf_id {0};
   int fission_nuclide;
   int surf_id {SURFACE_NONE};
   ParticleType particle;
@@ -434,11 +433,13 @@ class ParticleData : public GeometryState {
   double E_;
   double E_last_;
   double E_parent_; //!< energy of parent neutron in eV
+  double E_pprod_;
   int g_ {0};
   int g_last_;
 
   // Other birth data
   int fission_nuclide_;  //!< this particle was born as a result of this nuclide fissioning
+  int pprod_nuclide_;
   // a double for fission cross section at birth? if so, I need to also add it to the bank...
 
   double wgt_ {1.0};
@@ -480,6 +481,9 @@ class ParticleData : public GeometryState {
   vector<double> flux_derivs_;
 
   std::vector<std::vector<double>> cumulative_sensitivities_;  // for sensitivities for this particle
+  std::vector<std::vector<double>> cum_sens_;
+  std::vector<std::vector<double>> pprod_sens_;
+  double neutron_flux_last_ = 0.0;
 
   vector<FilterMatch> filter_matches_;
 
@@ -544,6 +548,8 @@ class ParticleData : public GeometryState {
   const double& E_last() const { return E_last_; }
   double& E_parent() { return E_parent_; }              // for sensitivity analysis
   const double& E_parent() const { return E_parent_; }  // for SA
+  double& E_pprod() { return E_pprod_; }
+  const double& E_pprod() const { return E_pprod_; }
   int& g() { return g_; }
   const int& g() const { return g_; }
   int& g_last() { return g_last_; }
@@ -645,8 +651,17 @@ class ParticleData : public GeometryState {
   // Used in sensitivity analysis
   std::vector<double>& cumulative_sensitivities(int i) { return cumulative_sensitivities_[i]; }
   const std::vector<double>& cumulative_sensitivities(int i) const { return cumulative_sensitivities_[i]; }
+  std::vector<double>& cum_sens(int i) { return cum_sens_[i]; }
+  const std::vector<double>& cum_sens(int i) const { return cum_sens_[i]; }
+  std::vector<double>& pprod_sens(int i) { return pprod_sens_[i]; }
+  const std::vector<double>& pprod_sens(int i) const { return pprod_sens_[i]; }
+  double neutron_flux_last() const { return neutron_flux_last_; }
+
   int& fission_nuclide() { return fission_nuclide_; }
   const int& fission_nuclide() const { return fission_nuclide_; }
+  
+  int& pprod_nuclide() { return pprod_nuclide_; }
+  const int& pprod_nuclide() const { return pprod_nuclide_; }
 
   // Matches of tallies
   decltype(filter_matches_)& filter_matches() { return filter_matches_; }
@@ -732,6 +747,26 @@ class ParticleData : public GeometryState {
     cumulative_sensitivities_[indx].resize(newSize, 0.0);
   }
 
+  void resize_init_cum_sens(int newSize)
+  {
+    cum_sens_.resize(newSize, {0.0});
+  }
+  
+  void resize_init_cum_sens_vec(int indx, int newSize)
+  {
+    cum_sens_[indx].resize(newSize, 0.0);
+  }
+  
+  void resize_init_pprod_sens(int newSize)
+  {
+    pprod_sens_.resize(newSize, {0.0});
+  }
+  
+  void resize_init_pprod_sens_vec(int indx, int newSize)
+  {
+    pprod_sens_[indx].resize(newSize, 0.0);
+  }
+  
   void resize_alloc_filter_matches(int newSize)
   {
     filter_matches_.resize(newSize);
@@ -749,6 +784,25 @@ class ParticleData : public GeometryState {
       std::fill(it.begin(), it.end(), 0.0);
     }
   }
+  
+  void initialize_cum_sens()
+  {
+    for (auto& it : cum_sens_){
+      std::fill(it.begin(), it.end(), 0.0);
+    }
+  }
+  
+  void initialize_pprod_sens()
+  {
+    for (auto& it : pprod_sens_){
+      std::fill(it.begin(), it.end(), 0.0);
+    }
+  }
+  
+  void set_neutron_flux_last(double v) 
+  { 
+    neutron_flux_last_ = v;
+  }
 };
 
 } // namespace openmc

include/openmc/simulation.h

@@ -75,6 +75,8 @@ void initialize_history(Particle& p, int64_t index_source);
 //! Helper function for initialize_history() that is called independently elsewhere
 void initialize_history_partial(Particle& p);
 
+void initialize_history_in_part(Particle& p);
+
 //! Finalize a batch
 //!
 //! Handles synchronization and accumulation of tallies, calculation of Shannon

include/openmc/tallies/sensitivity.h

@@ -74,7 +74,7 @@ struct TallySensitivity {
   int sens_reaction;    //!< Need something to specify reaction, use ReactionType?
   int sens_material;    //!< Material the sensitivity is calculated in
   int sens_cell;        //!< Cell the sensitivity is calculated in
-  std::vector<double> energy_bins_; //!< Energy bins on which to discretize the cross section
+  std::vector<double> energy_bins_; //!< Energy bins on which to tally the sensitivity
   int n_bins_; //!< something to indicate the size of the vector
 
   TallySensitivity() {}
@@ -112,7 +112,7 @@ void score_collision_sensitivity(Particle& p);
 //! \param distance The distance in [cm] traveled by the particle
 void score_track_sensitivity(Particle& p, double distance);
 
-void score_source_sensitivity(Particle& p);
+void score_pprod_sensitivity(Particle& p);
 
 } // namespace openmc
 

include/openmc/tallies/tally_scoring.h

@@ -99,7 +99,7 @@ void score_surface_tally(Particle& p, const vector<int>& tallies);
 
 //! Function to handle the special case when the tally is a sensitivity tally.
 void score_collision_sensitivity_tally(Particle& p, int i_tally, int start_index, int filter_index,
-  double filter_weight, int i_nuclide, double atom_density, double flux);
+  double filter_weight, int i_nuclide, double atom_density, double flux, double neutron_flux);
 
 //! Score the pulse-height tally
 //! This is triggered at the end of every particle history

src/particle.cpp

@@ -96,7 +96,6 @@ bool Particle::create_secondary(
   }
   bank.time = time();
   bank_second_E() += bank.E;
-  // if (!model::active_tallies.empty()) score_source_sensitivity(*this);
   return true;
 }
 
@@ -130,7 +129,7 @@ void Particle::from_source(const SourceSite* src)
   fission() = false;
   zero_flux_derivs();
 
-  initialize_cumulative_sensitivities();
+  initialize_cumulative_sensitivities();  
   lifetime() = 0.0;
 
   // Copy attributes from source bank site
@@ -163,6 +162,12 @@ void Particle::from_source(const SourceSite* src)
     int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
     surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
   }
+  
+  if (src->particle == ParticleType::neutron) {
+    initialize_cum_sens();
+    initialize_pprod_sens();
+  }
+  
 }
 
 void Particle::event_calculate_xs()
@@ -418,7 +423,7 @@ void Particle::event_collide()
   if (!model::active_tallies.empty()) {
     score_collision_derivative(*this);
     score_collision_sensitivity(*this);  // Score cumulative sensitivity for sensitivity tallies.
-    // score_source_sensitivity(*this);
+    // score_pprod_sensitivity(*this);
   }
 
 #ifdef DAGMC
@@ -476,7 +481,7 @@ void Particle::event_revive_from_secondary()
         n_coord_last() = n_coord();
       }
       pht_secondary_particles();
-    }
+    }    
 
     // Enter new particle in particle track file
     if (write_track())

src/simulation.cpp

@@ -641,6 +641,7 @@ void initialize_history(Particle& p, int64_t index_source)
     add_particle_track(p);
 
   initialize_history_partial(p);
+  initialize_history_in_part(p);
 }
 
 void initialize_history_partial(Particle& p)
@@ -664,6 +665,33 @@ void initialize_history_partial(Particle& p)
   p.resize_alloc_filter_matches(model::tally_filters.size());  
 }
 
+void initialize_history_in_part(Particle& p)
+{
+  
+  if (p.type() == ParticleType::photon) return;
+  
+  // for n-p coupled simulation
+  // and sensitivity of photon responses (heating, dose) to neutron cross-sections
+  // initialized variable to store derivatives of neutron transport operator 
+  // needed for photon transport sensitivity to neutron data
+  if (!model::active_tallies.empty())
+  {
+       
+    p.resize_init_cum_sens(model::tally_sens.size());
+    for (int idx=0; idx< model::tally_sens.size();idx++){
+      p.resize_init_cum_sens_vec(idx, model::tally_sens[idx].n_bins_);
+    }
+    
+    p.resize_init_pprod_sens(model::tally_sens.size());
+    for (int idx=0; idx< model::tally_sens.size();idx++){
+      p.resize_init_pprod_sens_vec(idx, model::tally_sens[idx].n_bins_);
+    }
+  }
+
+  // Allocate space for tally filter matches
+  // p.resize_alloc_filter_matches(model::tally_filters.size());  
+}
+
 int overall_generation()
 {
   using namespace simulation;