Implement relativistic point detector tally

Author: Itay

.gitmodules

@@ -13,3 +13,6 @@
 [submodule "vendor/Catch2"]
 	path = vendor/Catch2
 	url = https://github.com/catchorg/Catch2.git
+[submodule "vendor/gsl-lite"]
+    path = vendor/gsl-lite
+    url = https://github.com/gsl-lite/gsl-lite.git

include/openmc/constants.h

@@ -290,7 +290,7 @@ enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE };
 
 enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 
-enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, POINT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 

include/openmc/distribution.h

@@ -23,7 +23,7 @@ class Distribution {
 public:
   virtual ~Distribution() = default;
   virtual double sample(uint64_t* seed) const = 0;
-
+  virtual double get_pdf(double x) const = 0;
   //! Return integral of distribution
   //! \return Integral of distribution
   virtual double integral() const { return 1.0; };
@@ -84,6 +84,10 @@ class Discrete : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  //! Calculate the probability density function (PDF) at a given value
+  //! \param x The value at which to evaluate the PDF
+  //! \return The value of the PDF at the given point
+  double get_pdf(double x) const;
 
   double integral() const override { return di_.integral(); };
 
@@ -111,6 +115,7 @@ class Uniform : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -135,6 +140,7 @@ class PowerLaw : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
@@ -160,6 +166,7 @@ class Maxwell : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double theta() const { return theta_; }
 
@@ -180,6 +187,7 @@ class Watt : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double a() const { return a_; }
   double b() const { return b_; }
@@ -204,6 +212,7 @@ class Normal : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double mean_value() const { return mean_value_; }
   double std_dev() const { return std_dev_; }
@@ -227,8 +236,10 @@ class Tabular : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
-  // properties
+  // double get_pdf_value(double x,uint64_t* seed) const;
+  //  properties
   vector<double>& x() { return x_; }
   const vector<double>& x() const { return x_; }
   const vector<double>& p() const { return p_; }
@@ -263,6 +274,7 @@ class Equiprobable : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   const vector<double>& x() const { return x_; }
 
@@ -282,6 +294,7 @@ class Mixture : public Distribution {
   //! \param seed Pseudorandom number seed pointer
   //! \return Sampled value
   double sample(uint64_t* seed) const override;
+  double get_pdf(double x) const;
 
   double integral() const override { return integral_; }
 

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(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_;

include/openmc/hdf5_interface.h

@@ -114,7 +114,7 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name,
 void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
   const char* name, const long long* buffer, bool indep);
 void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
-  const char* name, char const* buffer, bool indep);
+  const char* name, const char* buffer, bool indep);
 void write_tally_results(
   hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results);
 } // extern "C"

include/openmc/particle.h

@@ -34,13 +34,14 @@ class Particle : public ParticleData {
 
   Particle() = default;
 
-  //==========================================================================
-  // Methods
-
   double speed() const;
 
   //! create a secondary particle
   //
+  //
+  double getMass() const;
+  //
+  //
   //! stores the current phase space attributes of the particle in the
   //! secondary bank and increments the number of sites in the secondary bank.
   //! \param wgt Weight of the secondary particle
@@ -63,7 +64,9 @@ class Particle : public ParticleData {
   //! simply as a secondary particle.
   //! \param src Source site data
   void from_source(const SourceSite* src);
-
+  void initialize_ghost_particle(Particle& p, Direction u, double E);
+  void initialize_ghost_particle_from_source(
+    const SourceSite* src, Direction u_new);
   // Coarse-grained particle events
   void event_calculate_xs();
   void event_advance();

include/openmc/particle_data.h

@@ -54,6 +54,8 @@ struct SourceSite {
   int parent_nuclide {-1};
   int64_t parent_id;
   int64_t progeny_id;
+  int source_index;
+  bool ext {true}; //!< Is this source site external?
 };
 
 //! State of a particle used for particle track files
@@ -488,7 +490,9 @@ class ParticleData : public GeometryState {
   int event_nuclide_;
   int event_mt_;
   int delayed_group_ {0};
+  int event_index_mt_;
   int parent_nuclide_ {-1};
+  Direction v_t_;
 
   int n_bank_ {0};
   double bank_second_E_ {0.0};
@@ -622,7 +626,10 @@ class ParticleData : public GeometryState {
   const int& event_nuclide() const { return event_nuclide_; }
   int& event_mt() { return event_mt_; } // MT number of collision
   const int& event_mt() const { return event_mt_; }
+  int& event_index_mt() { return event_index_mt_; }
   int& delayed_group() { return delayed_group_; } // delayed group
+  Position& v_t() { return v_t_; }
+  const Position& v_t() const { return v_t_; }
   const int& parent_nuclide() const { return parent_nuclide_; }
   int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide
 

include/openmc/physics.h

@@ -90,6 +90,11 @@ Direction sample_cxs_target_velocity(
 void sample_fission_neutron(
   int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p);
 
+void score_fission_neutron(int i_tally, int i_nuclide, const Reaction& rx,
+  SourceSite* site, Particle& p, std::vector<double>& mu_cm,
+  std::vector<double>& Js, std::vector<Particle>& ghost_particles,
+  std::vector<double>& pdfs_lab);
+
 //! handles all reactions with a single secondary neutron (other than fission),
 //! i.e. level scattering, (n,np), (n,na), etc.
 void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p);

include/openmc/reaction_product.h

@@ -11,8 +11,9 @@
 #include "openmc/endf.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/tally.h"
 #include "openmc/vector.h" // for vector
-
 namespace openmc {
 
 //==============================================================================
@@ -48,7 +49,10 @@ class ReactionProduct {
   //! \param[out] mu Outgoing cosine with respect to current direction
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
-
+  void get_pdf(int i_tally, double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

include/openmc/secondary_correlated.h

@@ -10,6 +10,7 @@
 #include "openmc/angle_energy.h"
 #include "openmc/distribution.h"
 #include "openmc/endf.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -40,6 +41,10 @@ class CorrelatedAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed,
+    Particle& p, std::vector<double>& mu_cm, std::vector<double>& Js,
+    std::vector<Particle>& ghost_particles,
+    std::vector<double>& pdfs_lab) const;
 
   // energy property
   vector<double>& energy() { return energy_; }