Point Detector Tally

include/openmc/constants.h

@@ -281,7 +281,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

@@ -22,7 +22,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; };
@@ -83,6 +83,7 @@ class Discrete : 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 di_.integral(); };
 
@@ -110,6 +111,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_; }
@@ -134,6 +136,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_); }
@@ -159,6 +162,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_; }
 
@@ -179,6 +183,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_; }
@@ -203,6 +208,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_; }
@@ -226,7 +232,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;
 
+
+  //double get_pdf_value(double x,uint64_t* seed) const;
   // properties
   vector<double>& x() { return x_; }
   const vector<double>& x() const { return x_; }
@@ -262,6 +271,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_; }
 
@@ -281,6 +291,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/particle.h

@@ -18,6 +18,7 @@
 
 namespace openmc {
 
+
 // Forward declare the Surface class for use in Particle::cross_vacuum_bc, etc.
 class Surface;
 
@@ -35,9 +36,6 @@ class Particle : public ParticleData {
 
   Particle() = default;
 
-  //==========================================================================
-  // Methods
-
   double speed() const;
 
   //! moves the particle by the distance length to its next location
@@ -46,6 +44,10 @@ class Particle : public ParticleData {
 
   //! 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
@@ -61,7 +63,8 @@ class Particle : public ParticleData {
   //! simply as a secondary particle.
   //! \param src Source site data
   void from_source(const SourceSite* src);
-
+  void initilze_ghost_particle(Particle& p,Direction u, double E);
+  void initilze_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

@@ -51,6 +51,7 @@ struct SourceSite {
   double wgt {1.0};
   int delayed_group {0};
   int surf_id {0};
+  bool ext {false};
   ParticleType particle;
   int64_t parent_id;
   int64_t progeny_id;
@@ -234,6 +235,7 @@ class ParticleData {
   //----------------------------------------------------------------------------
   // Constructors
   ParticleData();
+  void reset_cords();
 
 private:
   //==========================================================================
@@ -281,7 +283,9 @@ class ParticleData {
   TallyEvent event_;      //!< scatter, absorption
   int event_nuclide_;     //!< index in nuclides array
   int event_mt_;          //!< reaction MT
+  int event_index_mt_;
   int delayed_group_ {0}; //!< delayed group
+  Direction v_t_;
 
   // Post-collision physical data
   int n_bank_ {0};        //!< number of fission sites banked
@@ -420,7 +424,10 @@ class ParticleData {
   int& event_nuclide() { return event_nuclide_; }
   const int& event_nuclide() const { return event_nuclide_; }
   int& event_mt() { return event_mt_; }
+  int& event_index_mt() { return event_index_mt_; }
   int& delayed_group() { return delayed_group_; }
+  Position& v_t() { return v_t_; }
+  const Position& v_t() const { return v_t_; }
 
   int& n_bank() { return n_bank_; }
   int& n_bank_second() { return n_bank_second_; }

include/openmc/physics.h

@@ -90,6 +90,10 @@ 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,7 +11,8 @@
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
 #include "openmc/vector.h" // for vector
-
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/tally.h"
 namespace openmc {
 
 //==============================================================================
@@ -43,7 +44,7 @@ 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

@@ -11,6 +11,7 @@
 #include "openmc/distribution.h"
 #include "openmc/endf.h"
 #include "openmc/vector.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -40,6 +41,8 @@ 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_; }

include/openmc/secondary_kalbach.h

@@ -11,6 +11,7 @@
 #include "openmc/constants.h"
 #include "openmc/endf.h"
 #include "openmc/vector.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -31,6 +32,7 @@ class KalbachMann : 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;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

include/openmc/secondary_nbody.h

@@ -7,6 +7,7 @@
 #include "hdf5.h"
 
 #include "openmc/angle_energy.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -27,6 +28,9 @@ class NBodyPhaseSpace : 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;
+
 
 private:
   int n_bodies_;      //!< Number of particles distributed

include/openmc/secondary_thermal.h

@@ -32,7 +32,8 @@ class CoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
-
+ double get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const;
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
 };
@@ -55,6 +56,8 @@ class IncoherentElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   double debye_waller_;
@@ -80,6 +83,8 @@ class IncoherentElasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(
+    double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +111,8 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(
+    double E_in, double& E_out, double& mu, uint64_t* seed , int l = -1) const;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +141,8 @@ class IncoherentInelasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double get_pdf(
+  double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
 private:
   //! Secondary energy/angle distribution
@@ -178,6 +187,8 @@ class MixedElasticAE : public AngleEnergy {
   const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
 };
 
+
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H

include/openmc/secondary_uncorrelated.h

@@ -11,6 +11,7 @@
 #include "openmc/distribution_energy.h"
 #include "openmc/memory.h"
 #include "openmc/vector.h"
+#include "openmc/particle.h"
 
 namespace openmc {
 
@@ -31,6 +32,8 @@ class UncorrelatedAngleEnergy : 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;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

include/openmc/tallies/tally.h

@@ -55,6 +55,10 @@ class Tally {
 
   void set_nuclides(const vector<std::string>& nuclides);
 
+  void set_positions(pugi::xml_node node);
+
+  void set_positions(const vector<std::string>& positions);
+
   const xt::xtensor<double, 3>& results() const { return results_; }
 
   //! returns vector of indices corresponding to the tally this is called on
@@ -151,6 +155,9 @@ class Tally {
 
   vector<int> scores_; //!< Filter integrands (e.g. flux, fission)
 
+  //! Index of each pos to be tallied.
+  vector<std::string> positions_;
+
   //! Index of each nuclide to be tallied.  -1 indicates total material.
   vector<int> nuclides_ {-1};
 
@@ -206,6 +213,7 @@ extern vector<int> active_tallies;
 extern vector<int> active_analog_tallies;
 extern vector<int> active_tracklength_tallies;
 extern vector<int> active_collision_tallies;
+extern vector<int> active_point_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
 extern vector<int> active_pulse_height_tallies;