Testing point detectors

CMakeLists.txt

@@ -425,6 +425,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_musurface.cpp
   src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
+  src/tallies/filter_point.cpp  
   src/tallies/filter_polar.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp

include/openmc/angle_energy.h

@@ -16,6 +16,8 @@ class AngleEnergy {
 public:
   virtual void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
+  virtual double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 

include/openmc/chain.h

@@ -70,6 +70,8 @@ class DecayPhotonAngleEnergy : public AngleEnergy {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const Distribution* photon_energy_;

include/openmc/constants.h

@@ -294,9 +294,15 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
-enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
+enum class TallyType {
+  VOLUME,
+  MESH_SURFACE,
+  SURFACE,
+  PULSE_HEIGHT,
+  NEXT_EVENT
+};
 
-enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
+enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, POINT };
 
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 

include/openmc/distribution_angle.h

@@ -25,6 +25,7 @@ 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 evaluate(double E, double mu) const;
 
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty

include/openmc/distribution_multi.h

@@ -6,6 +6,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/distribution.h"
+#include "openmc/error.h"
 #include "openmc/position.h"
 
 namespace openmc {
@@ -29,6 +30,11 @@ class UnitSphereDistribution {
   //! \return (sampled Direction, sample weight)
   virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
 
+  virtual double evaluate(Direction u) const
+  {
+    fatal_error("evaluate not available for this UnitSphereDistribution type");
+  }
+
   Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
 };
 
@@ -52,6 +58,8 @@ class PolarAzimuthal : public UnitSphereDistribution {
   //! \return (sampled Direction, value of the PDF at this Direction)
   std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
 
+  double evaluate(Direction u) const override;
+
   // Observing pointers
   Distribution* mu() const { return mu_.get(); }
   Distribution* phi() const { return phi_.get(); }
@@ -87,6 +95,8 @@ class Isotropic : public UnitSphereDistribution {
   //! \return (sampled direction, sample weight)
   std::pair<Direction, double> sample(uint64_t* seed) const override;
 
+  double evaluate(Direction u) const override;
+
   // Set or get bias distribution
   void set_bias(std::unique_ptr<PolarAzimuthal> bias)
   {

include/openmc/math_functions.h

@@ -201,6 +201,12 @@ std::complex<double> faddeeva(std::complex<double> z);
 //! \return Derivative of Faddeeva function evaluated at z
 std::complex<double> w_derivative(std::complex<double> z, int order);
 
+//! Evaluate relative exponential function
+//!
+//! \param x Real argument
+//! \return (exp(x)-1)/x without loss of precision near 0
+double exprel(double x);
+
 //! Helper function to get index and interpolation function on an incident
 //! energy grid
 //!

include/openmc/particle.h

@@ -37,6 +37,8 @@ class Particle : public ParticleData {
   //==========================================================================
   // Methods
 
+  double mass() const;
+
   double speed() const;
 
   //! create a secondary particle
@@ -64,10 +66,12 @@ class Particle : public ParticleData {
   //! \param src Source site data
   void from_source(const SourceSite* src);
 
+  void initialize_pseudoparticle(Particle& p, Direction u, double E);
+
   // Coarse-grained particle events
-  void event_calculate_xs();
+  void event_calculate_xs(bool is_pseudo = false);
   void event_advance();
-  void event_cross_surface();
+  void event_cross_surface(bool is_pseudo = false);
   void event_collide();
   void event_revive_from_secondary();
   void event_death();

include/openmc/physics.h

@@ -85,7 +85,7 @@ void sample_fission_neutron(
 
 //! 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);
+void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p);
 
 void sample_secondary_photons(Particle& p, int i_nuclide);
 

include/openmc/position.h

@@ -204,6 +204,17 @@ inline Position operator/(double a, Position b)
   return b /= a;
 }
 
+inline bool operator<(Position a, Position b)
+{
+  if (a[0] != b[0])
+    return (a[0] < b[0]);
+  if (a[1] != b[1])
+    return (a[1] < b[1]);
+  if (a[2] != b[2])
+    return (a[2] < b[2]);
+  return false;
+}
+
 inline Position Position::reflect(Position n) const
 {
   const double projection = n.dot(*this);

include/openmc/random_lcg.h

@@ -9,11 +9,12 @@ namespace openmc {
 // Module constants.
 //==============================================================================
 
-constexpr int N_STREAMS {4};
+constexpr int N_STREAMS {5};
 constexpr int STREAM_TRACKING {0};
 constexpr int STREAM_SOURCE {1};
 constexpr int STREAM_URR_PTABLE {2};
 constexpr int STREAM_VOLUME {3};
+constexpr int STREAM_NEXT_EVENT {4};
 constexpr int64_t DEFAULT_SEED {1};
 constexpr uint64_t DEFAULT_STRIDE {152917ULL};
 

include/openmc/reaction_product.h

@@ -49,6 +49,11 @@ class ReactionProduct {
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) 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

@@ -40,6 +40,9 @@ 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;
+  Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // energy property
   vector<double>& energy() { return energy_; }

include/openmc/secondary_kalbach.h

@@ -31,6 +31,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 sample_params(double E_in, double& E_out, double& km_a, double& km_r,
+    uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Outgoing energy/angle at a single incoming energy

include/openmc/secondary_nbody.h

@@ -27,6 +27,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;
+  double sample_energy(double E_in, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   int n_bodies_;      //!< Number of particles distributed

include/openmc/secondary_thermal.h

@@ -6,6 +6,7 @@
 
 #include "openmc/angle_energy.h"
 #include "openmc/endf.h"
+#include "openmc/search.h"
 #include "openmc/secondary_correlated.h"
 #include "openmc/vector.h"
 
@@ -32,6 +33,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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@@ -55,6 +58,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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   double debye_waller_;
@@ -80,6 +85,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 sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
@@ -106,6 +113,9 @@ 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;
+  void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -134,6 +144,10 @@ 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;
+  void sample_params(double E_in, double& E_out, double& f, int& l, int& j,
+    uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   //! Secondary energy/angle distribution
@@ -169,6 +183,9 @@ class MixedElasticAE : public AngleEnergy {
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
@@ -178,6 +195,52 @@ class MixedElasticAE : public AngleEnergy {
   const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
 };
 
+struct DoubleVector {
+  double data;
+  const double& operator[](size_t index) const { return data; }
+};
+
+template<typename T>
+double get_pdf_discrete(const vector<double>& mu, const T& w, double mu_0)
+{
+  // Make sure mu is in range [-1,1]
+  if (std::abs(mu_0) > 1.0)
+    mu_0 = std::copysign(1.0, mu_0);
+  double a0;
+  double a1;
+  double b0;
+  double b1;
+  int32_t ai = -1;
+  int32_t bi = -1;
+  if (mu_0 > mu[0]) {
+    ai = lower_bound_index(mu.begin(), mu.end(), mu_0);
+    a0 = mu[ai];
+    a1 = (ai > 1) ? mu[ai - 1] : -1.0;
+  } else {
+    a0 = -1.0;
+    a1 = -1.0;
+  }
+  if (mu_0 < mu[mu.size() - 1]) {
+    bi = upper_bound_index(mu.begin(), mu.end(), mu_0);
+    b0 = mu[bi];
+    b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;
+  } else {
+    b0 = 1.0;
+    b1 = 1.0;
+  }
+
+  //  Calculate Delta_a and Delta_b
+  double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);
+  double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);
+
+  if (mu_0 < a0 + delta_a)
+    return w[ai] / (2.0 * delta_a);
+  else if (mu_0 + delta_b < b0)
+    return w[bi] / (2.0 * delta_b);
+  else
+    return 0.0;
+}
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H

include/openmc/secondary_uncorrelated.h

@@ -31,6 +31,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;
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
   // Accessors
   AngleDistribution& angle() { return angle_; }

include/openmc/simulation.h

@@ -36,12 +36,14 @@ extern "C" double
 extern double log_spacing; //!< lethargy spacing for energy grid searches
 extern "C" int n_lost_particles;   //!< cumulative number of lost particles
 extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
-extern "C" int restart_batch;      //!< batch at which a restart job resumed
-extern "C" bool satisfy_triggers;  //!< have tally triggers been satisfied?
-extern int ssw_current_file;       //!< current surface source file
-extern "C" int total_gen;          //!< total number of generations simulated
-extern double total_weight;        //!< Total source weight in a batch
-extern int64_t work_per_rank;      //!< number of particles per MPI rank
+extern bool
+  nonvacuum_boundary_present;     //!< Does the geometry contain non-vacuum b.c.
+extern "C" int restart_batch;     //!< batch at which a restart job resumed
+extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
+extern int ssw_current_file;      //!< current surface source file
+extern "C" int total_gen;         //!< total number of generations simulated
+extern double total_weight;       //!< Total source weight in a batch
+extern int64_t work_per_rank;     //!< number of particles per MPI rank
 
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
@@ -92,6 +94,9 @@ void broadcast_results();
 
 void free_memory_simulation();
 
+//! Simulate a single pseudoparticle history
+double transport_pseudoparticle(Particle& p, double total_distance);
+
 //! Simulate a single particle history (and all generated secondary particles,
 //!  if enabled), from birth to death
 void transport_history_based_single_particle(Particle& p);

include/openmc/tallies/filter.h

@@ -39,6 +39,7 @@ enum class FilterType {
   MUSURFACE,
   PARENT_NUCLIDE,
   PARTICLE,
+  POINT,
   POLAR,
   SPHERICAL_HARMONICS,
   SPATIAL_LEGENDRE,