created get pdf value and didnt use it

openmc-dev · Oct 25, 2022 · 8269657 · Itay-max · Oct 26, 2022 · Itay-max
1 parent 21efc8d
commit 8269657
Show file tree

Hide file tree

Showing 6 changed files with 93 additions and 13 deletions.
diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
@@ -215,7 +215,7 @@ class Tabular : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const;
-
+  double get_pdf_value(double x,uint64_t* seed) const;
   // x property
   vector<double>& x() { return x_; }
   const vector<double>& x() const { return x_; }

diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h
@@ -25,10 +25,11 @@ class AngleDistribution {
   //! \param[inout] seed pseudorandom number seed pointer
   //! \return Cosine of the angle in the range [-1,1]
   double sample(double E, uint64_t* seed) const;
-
+  double get_pdf_value(double E,double mu, uint64_t* seed) const;
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
   bool empty() const { return energy_.empty(); }
+  double get_energy(int num) {return energy_[num];}
 
 private:
   vector<double> energy_;

diff --git a/src/distribution.cpp b/src/distribution.cpp
@@ -14,6 +14,7 @@
 #include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
+#include "openmc/search.h"
 
 namespace openmc {
 
@@ -299,6 +300,42 @@ double Tabular::sample(uint64_t* seed) const
   }
 }
 
+double Tabular::get_pdf_value(double x,uint64_t* seed) const
+{
+  // get PDF value at x 
+
+  int i;
+  std::size_t n = x_.size();
+  if (x < x_[0]) {
+    return 0;
+  } else if (x > x_[n - 1]) {
+    return 0;
+  } else {
+    i = lower_bound_index(x_.begin(), x_.end(), x);
+  }
+
+  // Determine bounding PDF values
+  double x_i = x_[i];
+  double p_i = p_[i];
+
+  if (interp_ == Interpolation::histogram) {
+    // Histogram interpolation
+    return p_i;
+  } else {
+    // Linear-linear interpolation
+    double x_i1 = x_[i + 1];
+    double p_i1 = p_[i + 1];
+
+    double m = (p_i1 - p_i) / (x_i1 - x_i);
+    if (m == 0.0) {
+      return p_i;
+    } else {
+      return p_i + (x-x_i)*m;
+    }
+  }
+}
+
+
 //==============================================================================
 // Equiprobable implementation
 //==============================================================================

diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp
@@ -95,4 +95,37 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::get_pdf_value(double E,double mu, uint64_t* seed) const
+{
+  // Determine number of incoming energies
+  auto n = energy_.size();
+
+  // Find energy bin and calculate interpolation factor -- if the energy is
+  // outside the range of the tabulated energies, choose the first or last bins
+  int i;
+  double r;
+  if (E < energy_[0]) {
+    i = 0;
+    r = 0.0;
+  } else if (E > energy_[n - 1]) {
+    i = n - 2;
+    r = 1.0;
+  } else {
+    i = lower_bound_index(energy_.begin(), energy_.end(), E);
+    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
+  }
+
+  // Sample between the ith and (i+1)th bin
+  if (r > prn(seed))
+    ++i;
+
+  // Sample i-th distribution
+  double pdf = distribution_[i]->get_pdf_value(mu,seed);
+
+  // Make sure pdf > 0 and return
+  if (std::abs(pdf) < 0)
+    fmt::print("pdf = {} <1",pdf);
+  return pdf;
+}
+
 } // namespace openmc
diff --git a/src/particle.cpp b/src/particle.cpp
@@ -113,6 +113,12 @@ void Particle::from_source(const SourceSite* src)
   E_last() = E();
   time() = src->time;
   time_last() = src->time;
+
+  fmt::print("==============================particle created==============================\n");
+  //fmt::print("u = {0} , {1} , {2}\n",u().x,u().y,u().z);
+  //fmt::print("u_last = {0} , {1} , {2}\n",u_last().x,u_last().y,u_last().z);
+  //fmt::print("pos = {0} , {1} , {2}\n",r().x,r().y,r().z);
+
 }
 
 void Particle::event_calculate_xs()

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
@@ -2370,27 +2370,30 @@ void score_collision_tally(Particle& p)
   //fmt::print("collison happened\n");
   if(p.E_last()>2e5)
   {
-  fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last());
-  fmt::print("event mt = {}\n",p.event_mt());
-  fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last());
+  //fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last());
+  //fmt::print("event mt = {}\n",p.event_mt());
+  //fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last());
   //fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z);
-
-  fmt::print("distance = {}\n",distance);
-  fmt::print("xs = {}\n",p.macro_xs().total);
+  //fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z);
+  //fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z);
+  //fmt::print("distance = {}\n",distance);
+  //fmt::print("xs = {}\n",p.macro_xs().total);
 
   const auto& mat {model::materials[p.material()]};
   int n = mat->nuclide_.size();
   const auto& nuc {data::nuclides[0]};
   const auto& nuc2 {data::nuclides[1]};
-  fmt::print("number of nuclear = {}\n",n);
-  fmt::print("name of 1# nuclicde = {}\n",nuc->name_);
-  fmt::print("name of 2# nuclicde = {}\n",nuc2->name_);
+  //fmt::print("number of nuclear = {}\n",n);
+  //fmt::print("name of 1# nuclicde = {}\n",nuc->name_);
+  //fmt::print("name of 2# nuclicde = {}\n",nuc2->name_);
   int i_nuclide=0;
   const auto& rx {nuc->reactions_[0]};
+  //fmt::print("reaction MT = {}\n",rx->mt_);
   auto& d = rx->products_[0].distribution_[0];
   auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
-
-  fmt::print("name of 2# nuclicde = {}\n",d_.energy_[0]);
+  //fmt::print("d type = {}\n",typeid(d.get()).name());
+  //fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0));
+  //fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1));
   }
   //fmt::print("post energy = {}\n",p.E);
   //fmt::print(p.E);