Point Detector Tally

.clang-format

@@ -84,7 +84,6 @@ PenaltyBreakTemplateDeclaration: 10
 PenaltyExcessCharacter: 1000000
 PenaltyReturnTypeOnItsOwnLine: 200
 PointerAlignment: Left
-QualifierAlignment: Left
 ReflowComments:  true
 SortIncludes:    true
 SortUsingDeclarations: true

Dockerfile

@@ -182,8 +182,8 @@ FROM dependencies AS build
 
 ENV HOME=/root
 
-ARG openmc_branch=master
-ENV OPENMC_REPO='https://github.com/openmc-dev/openmc'
+ARG openmc_branch=deploy
+ENV OPENMC_REPO='https://github.com/itay-space/openmc.git'
 
 ARG compile_cores
 ARG build_dagmc

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

@@ -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,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(); };
 
@@ -111,6 +112,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 +137,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 +163,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 +184,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 +209,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 +233,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 +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_; }
 
@@ -282,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/ncrystal_interface.h

@@ -58,19 +58,10 @@ class NCrystalMat {
 
   //----------------------------------------------------------------------------
   // Trivial methods when compiling without NCRYSTAL
-  std::string cfg() const
-  {
-    return "";
-  }
-  double xs(const Particle& p) const
-  {
-    return -1.0;
-  }
+  std::string cfg() const { return ""; }
+  double xs(const Particle& p) const { return -1.0; }
   void scatter(Particle& p) const {}
-  operator bool() const
-  {
-    return false;
-  }
+  operator bool() const { return false; }
 #endif
 
 private:

include/openmc/particle.h

@@ -34,9 +34,6 @@ class Particle : public ParticleData {
 
   Particle() = default;
 
-  //==========================================================================
-  // Methods
-
   double speed() const;
 
   //! moves the particle by the distance length to its next location
@@ -45,6 +42,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
@@ -60,7 +61,9 @@ 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

@@ -48,9 +48,11 @@ 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;
+  int source_index;
 };
 
 //! State of a particle used for particle track files
@@ -391,41 +393,78 @@ class GeometryState {
 class ParticleData : public GeometryState {
 private:
   //==========================================================================
-  // Data members -- see public: below for descriptions
+  // Data members (accessor methods are below)
 
-  vector<NuclideMicroXS> neutron_xs_;
-  vector<ElementMicroXS> photon_xs_;
-  MacroXS macro_xs_;
-  CacheDataMG mg_xs_cache_;
+  // Cross section caches
+  vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
+  vector<ElementMicroXS> photon_xs_;  //!< Microscopic photon cross sections
+  MacroXS macro_xs_;                  //!< Macroscopic cross sections
+  CacheDataMG mg_xs_cache_;           //!< Multigroup XS cache
 
-  ParticleType type_ {ParticleType::neutron};
+  int64_t id_;                                //!< Unique ID
+  ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.)
 
-  double E_;
-  double E_last_;
-  int g_ {0};
-  int g_last_;
+  int n_coord_ {1};          //!< number of current coordinate levels
+  int cell_instance_;        //!< offset for distributed properties
+  vector<LocalCoord> coord_; //!< coordinates for all levels
 
-  double wgt_ {1.0};
-  double mu_;
-  double time_ {0.0};
-  double time_last_ {0.0};
-  double wgt_last_ {1.0};
+  // Particle coordinates before crossing a surface
+  int n_coord_last_ {1};  //!< number of current coordinates
+  vector<int> cell_last_; //!< coordinates for all levels
+
+  // Energy data
+  double E_;      //!< post-collision energy in eV
+  double E_last_; //!< pre-collision energy in eV
+  int g_ {0};     //!< post-collision energy group (MG only)
+  int g_last_;    //!< pre-collision energy group (MG only)
 
-  bool fission_ {false};
-  TallyEvent event_;
-  int event_nuclide_;
-  int event_mt_;
-  int delayed_group_ {0};
+  // Other physical data
+  double wgt_ {1.0};       //!< particle weight
+  double mu_;              //!< angle of scatter
+  double time_ {0.0};      //!< time in [s]
+  double time_last_ {0.0}; //!< previous time in [s]
+
+  // Other physical data
+  Position r_last_current_; //!< coordinates of the last collision or
+                            //!< reflective/periodic surface crossing for
+                            //!< current tallies
+  Position r_last_;         //!< previous coordinates
+  Direction u_last_;        //!< previous direction coordinates
+  double wgt_last_ {1.0};   //!< pre-collision particle weight
+
+  // What event took place
+  bool fission_ {false}; //!< did particle cause implicit fission
+  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
+  int n_bank_second_ {0}; //!< number of secondary particles banked
+  double wgt_bank_ {0.0}; //!< weight of fission sites banked
+  int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission
+                                           //!< sites banked
+
+  // Indices for various arrays
+  int surface_ {0};        //!< index for surface particle is on
+  int cell_born_ {-1};     //!< index for cell particle was born in
+  int material_ {-1};      //!< index for current material
+  int material_last_ {-1}; //!< index for last material
 
-  int n_bank_ {0};
-  int n_bank_second_ {0};
-  double wgt_bank_ {0.0};
-  int n_delayed_bank_[MAX_DELAYED_GROUPS];
+  // Boundary information
+  BoundaryInfo boundary_;
 
-  int cell_born_ {-1};
+  // Temperature of current cell
+  double sqrtkT_ {-1.0};     //!< sqrt(k_Boltzmann * temperature) in eV
+  double sqrtkT_last_ {0.0}; //!< last temperature
 
-  int n_collision_ {0};
+  // Statistical data
+  int n_collision_ {0}; //!< number of collisions
 
+  // Track output
   bool write_track_ {false};
 
   uint64_t seeds_[N_STREAMS];
@@ -528,8 +567,11 @@ class ParticleData : public GeometryState {
   bool& fission() { return fission_; }            // true if implicit fission
   int& event_nuclide() { return event_nuclide_; } // index of collision nuclide
   const int& event_nuclide() const { return event_nuclide_; }
-  int& event_mt() { return event_mt_; }           // MT number of collision
+  int& event_mt() { return event_mt_; } // MT number of collision
+  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_; }
 
   // Post-collision data
   int& n_bank() { return n_bank_; } // number of banked fission sites

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

@@ -10,8 +10,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 {
 
 //==============================================================================
@@ -43,7 +44,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_; }

include/openmc/secondary_kalbach.h

@@ -10,6 +10,7 @@
 #include "openmc/angle_energy.h"
 #include "openmc/constants.h"
 #include "openmc/endf.h"
+#include "openmc/particle.h"
 #include "openmc/vector.h"
 
 namespace openmc {
@@ -31,6 +32,10 @@ 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