crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge bfaf966173 into 61c8a59cff

Browse source
This commit is contained in:
GuySten 2026-07-20 21:24:01 +03:00 • committed by GitHub
commit 0a0a838b6b
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
44 changed files with 837 additions and 110 deletions
Download patch file Download diff file Expand all files Collapse all files

3
CMakeLists.txt
View file

@ -355,6 +355,7 @@ endif()
#===============================================================================
list(APPEND libopenmc_SOURCES
src/angle_energy.cpp
src/atomic_mass.cpp
src/bank.cpp
src/boundary_condition.cpp
@ -451,6 +452,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_particle_production.cpp
src/tallies/filter_polar.cpp
src/tallies/filter_reaction.cpp
@ -461,6 +463,7 @@ list(APPEND libopenmc_SOURCES
src/tallies/filter_universe.cpp
src/tallies/filter_weight.cpp
src/tallies/filter_zernike.cpp
src/tallies/next_event_scoring.cpp
src/tallies/tally.cpp
src/tallies/tally_scoring.cpp
src/tallies/trigger.cpp

14
include/openmc/angle_energy.h
View file

@ -27,12 +27,20 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
virtual double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
virtual double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const = 0;
virtual ~AngleEnergy() = default;
};
double get_jac_and_transform(
double E_in, double& mu, double& E_out, uint64_t* seed, double awr);
double get_jac_and_transform_impl(
double E_com, double& mu, double& E_out, uint64_t* seed, double awr);
} // namespace openmc
#endif // OPENMC_ANGLE_ENERGY_H

8
include/openmc/chain.h
View file

@ -76,9 +76,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
const Distribution* photon_energy_;

4
include/openmc/constants.h
View file

@ -296,9 +296,9 @@ 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, POINT };
enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, NEXT_EVENT };
enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };

2
include/openmc/particle_data.h
View file

@ -487,7 +487,7 @@ private:
* Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
* https://doi.org/10.1016/j.anucene.2017.11.032.
*/
class ParticleData : public GeometryState {
class ParticleData : virtual public GeometryState {
private:
//==========================================================================
// Data members -- see public: below for descriptions

2
include/openmc/physics.h
View file

@ -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);

11
include/openmc/position.h
View file

@ -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);

35
include/openmc/ray.h
View file

@ -1,14 +1,14 @@
#ifndef OPENMC_RAY_H
#define OPENMC_RAY_H
#include "openmc/particle_data.h"
#include "openmc/particle.h"
#include "openmc/position.h"
namespace openmc {
// Base class that implements ray tracing logic, not necessarily through
// defined regions of the geometry but also outside of it.
class Ray : public GeometryState {
class Ray : virtual public GeometryState {
public:
// Initialize from location and direction
@ -24,7 +24,7 @@ public:
* Traces the ray through the geometry, calling on_intersection
* at every surface boundary.
*/
void trace();
void trace(double max_distance = INFTY);
// Stops the ray and exits tracing when called from on_intersection
void stop() { stop_ = true; }
@ -32,6 +32,8 @@ public:
// Sets the dist_ variable
void compute_distance();
virtual void update_distance();
protected:
// Records how far the ray has traveled
double traversal_distance_ {0.0};
@ -46,5 +48,32 @@ private:
unsigned event_counter_ {0};
};
class ParticleRay : public Ray, public Particle {
public:
ParticleRay(
Position r, Direction u, ParticleType type_, double time_, double E_)
: Ray(r, u)
{
type() = type_;
time() = time_;
E() = E_;
E_last() = E_;
r_last() = r;
}
void on_intersection() override;
// Sets the dist_ variable
void update_distance() override;
const double& traversal_distance() const { return traversal_distance_; }
const double& traversal_mfp() const { return traversal_mfp_; }
protected:
// Records how much mean free paths the ray traveled
double traversal_mfp_ {0.0};
};
} // namespace openmc
#endif // OPENMC_RAY_H

8
include/openmc/reaction_product.h
View file

@ -60,9 +60,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const;
ParticleType particle_; //!< Particle type
EmissionMode emission_mode_; //!< Emission mode

8
include/openmc/secondary_correlated.h
View file

@ -53,9 +53,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
// energy property
vector<double>& energy() { return energy_; }

8
include/openmc/secondary_kalbach.h
View file

@ -46,9 +46,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
//! Outgoing energy/angle at a single incoming energy

8
include/openmc/secondary_nbody.h
View file

@ -39,9 +39,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
int n_bodies_; //!< Number of particles distributed

48
include/openmc/secondary_thermal.h
View file

@ -39,9 +39,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
@ -74,9 +76,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
double debye_waller_;
@ -108,9 +112,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
const vector<double>& energy_; //!< Energies at which cosines are tabulated
@ -149,9 +155,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
const vector<double>& energy_; //!< Incident energies
@ -196,9 +204,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
//! Secondary energy/angle distribution
@ -246,9 +256,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
private:
CoherentElasticAE coherent_dist_; //!< Coherent distribution

8
include/openmc/secondary_uncorrelated.h
View file

@ -37,9 +37,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const override;
// Accessors
AngleDistribution& angle() { return angle_; }

14
include/openmc/simulation.h
View file

@ -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;

1
include/openmc/tallies/filter.h
View file

@ -40,6 +40,7 @@ enum class FilterType {
PARENT_NUCLIDE,
PARTICLE,
PARTICLE_PRODUCTION,
POINT,
POLAR,
REACTION,
SPHERICAL_HARMONICS,

55
include/openmc/tallies/filter_point.h Normal file
View file

@ -0,0 +1,55 @@
#ifndef OPENMC_TALLIES_FILTER_POINT_H
#define OPENMC_TALLIES_FILTER_POINT_H
#include "openmc/position.h"
#include "openmc/span.h"
#include "openmc/tallies/filter.h"
#include "openmc/vector.h"
namespace openmc {
//==============================================================================
//! Bins tally by point detectors
//==============================================================================
class PointFilter : public Filter {
public:
//----------------------------------------------------------------------------
// Constructors, destructors
~PointFilter() = default;
//----------------------------------------------------------------------------
// Methods
std::string type_str() const override { return "point"; }
FilterType type() const override { return FilterType::POINT; }
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
//----------------------------------------------------------------------------
// Accessors
const vector<std::pair<Position, double>>& detectors() const
{
return detectors_;
}
void set_detectors(span<std::pair<Position, double>> detectors);
private:
//----------------------------------------------------------------------------
// Data members
vector<std::pair<Position, double>> detectors_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_PARTICLE_H

60
include/openmc/tallies/next_event_scoring.h Normal file
View file

@ -0,0 +1,60 @@
#ifndef OPENMC_TALLIES_NEXT_EVENT_SCORING_H
#define OPENMC_TALLIES_NEXT_EVENT_SCORING_H
#include "openmc/nuclide.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/ray.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/thermal.h"
namespace openmc {
void score_point_tally_elastic(
Particle& p, int i_nuclide, const Reaction& rx, int i_product, Direction v_t);
void score_point_tally_inelastic(
Particle& p, int i_nuclide, const Reaction& rx, int i_product, double yield);
void score_point_tally_fission(
Particle& p, int i_nuclide, const Reaction& rx, int i_product);
void score_point_tally_sab(Particle& p, int i_nuclide, const ThermalData& sab,
const NuclideMicroXS& micro);
void score_point_tally_source(SourceSite& site, int source_index);
template<typename PDF>
void score_point_tally_impl(
const Position r, const ParticleType type, const double time, PDF pdffunc)
{
for (auto& det : model::active_point_detectors) {
auto u = (det - r);
double total_distance = u.norm();
u /= total_distance;
double E;
double pdf = pdffunc(u, E);
if (pdf == 0.0)
continue;
auto p = ParticleRay(r, u, type, time, E);
p.Ray::trace(total_distance);
double distance = p.traversal_distance();
if (distance < total_distance)
continue;
double mfp = p.traversal_mfp();
double attenuation = std::exp(-mfp);
// Save the attenuation for point filter handling
p.wgt_last() = p.wgt();
p.wgt() *= attenuation;
double flux = p.wgt_last() * pdf;
score_tracklength_tally_general(p, flux, model::active_point_tallies);
}
}
} // namespace openmc
#endif // OPENMC_TALLIES_NEXT_EVENT_SCORING_H

4
include/openmc/tallies/tally.h
View file

@ -3,6 +3,7 @@
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/position.h"
#include "openmc/span.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/trigger.h"
@ -11,6 +12,7 @@
#include "openmc/tensor.h"
#include "pugixml.hpp"
#include <set>
#include <string>
#include <unordered_map>
@ -209,11 +211,13 @@ extern vector<int> active_analog_tallies;
extern vector<int> active_tracklength_tallies;
extern vector<int> active_timed_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;
extern vector<int32_t> pulse_height_cells;
extern vector<double> time_grid;
extern std::set<Position> active_point_detectors;
} // namespace model

6
include/openmc/tallies/tally_scoring.h
View file

@ -1,9 +1,12 @@
#ifndef OPENMC_TALLIES_TALLY_SCORING_H
#define OPENMC_TALLIES_TALLY_SCORING_H
#include "openmc/nuclide.h"
#include "openmc/particle.h"
#include "openmc/ray.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/thermal.h"
namespace openmc {
@ -81,6 +84,9 @@ void score_analog_tally_ce(Particle& p);
//! \param p The particle being tracked
void score_analog_tally_mg(Particle& p);
void score_tracklength_tally_general(
Particle& p, double flux, const vector<int>& tallies);
//! Score tallies using a tracklength estimate of the flux.
//
//! This is triggered at every event (surface crossing, lattice crossing, or

6
include/openmc/thermal.h
View file

@ -67,9 +67,11 @@ public:
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
//! \param[in] is_com Is scattering cosine is given in center of mass
//! coordinates \param[in] awr Weight of nucleus in neutron masses \return
//! Probability density for the scattering cosine
double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
double mu, double& E_out, uint64_t* seed) const;
double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const;
private:
struct Reaction {

84
openmc/filter.py
View file

@ -26,7 +26,7 @@ _FILTER_TYPES = (
'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell',
'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre',
'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle',
'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide',
'particleproduction', 'point', 'cellinstance', 'collision', 'time', 'parentnuclide',
'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
)
@ -793,6 +793,88 @@ class ParticleFilter(Filter):
return cls(bins, filter_id=filter_id)
class PointFilter(Filter):
"""Bins tally events based on point detectors.
Parameters
----------
bins : sequence of tuple[tuple[Real, Real, Real], Real]
Point detectors positions and exclusion radii.
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : sequence of tuple[tuple[Real, Real, Real], Real]
Point detectors positions and exclusion radii.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
"""
__hash__ = Filter.__hash__
def __eq__(self, other):
if type(self) is not type(other):
return False
elif len(self.bins) != len(other.bins):
return False
else:
return all(b1==b2 for b1,b2 in zip(self.bins,other.bins))
@Filter.bins.setter
def bins(self, bins):
cv.check_type('bins', bins, Sequence, tuple)
for i, item in enumerate(bins):
cv.check_type(f'bins[{i}]', item, tuple)
cv.check_length(f'bins[{i}]', item, 2, 2)
cv.check_type(f'bins[{i}][0]', item[0], tuple, Real)
cv.check_length(f'bins[{i}][0]', item[0], 3, 3)
cv.check_type(f'bins[{i}][1]', item[1], Real)
self._bins = bins
@classmethod
def from_hdf5(cls, group, **kwargs):
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
flat = group['bins'][()]
# Reconstruct tuple structure: every 4 values = (x, y, z, r0)
bins = []
for i in range(0, len(flat), 4):
pos = (float(flat[i]), float(flat[i+1]), float(flat[i+2]))
r0 = float(flat[i+3])
bins.append((pos, r0))
out = cls(bins, filter_id=filter_id)
out._num_bins = group['n_bins'][()]
return out
def get_pandas_dataframe(self, data_size, stride, **kwargs):
import pandas as pd
labels = [f"({p[0]}, {p[1]}, {p[2]}) R0={r}" for (p, r) in self.bins]
filter_bins = np.repeat(labels, stride)
tile_factor = data_size // len(filter_bins)
filter_bins = np.tile(filter_bins, tile_factor)
return pd.DataFrame({self.short_name.lower(): filter_bins})
def to_xml_element(self):
"""Return XML Element representing the Filter.
Returns
-------
element : lxml.etree._Element
XML element containing filter data
"""
element = ET.Element('filter')
element.set('id', str(self.id))
element.set('type', self.short_name.lower())
subelement = ET.SubElement(element, 'bins')
subelement.text = ' '.join(str(b) for item in self.bins
for b in list(item[0])+[item[1]])
return element
class ParentNuclideFilter(ParticleFilter):
"""Bins tally events based on the parent nuclide

4
openmc/lib/filter.py
View file

@ -22,7 +22,7 @@ __all__ = [
'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter',
'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter',
'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PointFilter', 'PolarFilter',
'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter',
'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter',
'ZernikeRadialFilter', 'filters'
@ -607,6 +607,8 @@ class ParticleFilter(Filter):
self._index, particle_i.ctypes.data_as(POINTER(c_int32)))
return [ParticleType(i) for i in particle_i]
class PointFilter(Filter):
filter_type = 'point'
class ParticleProductionFilter(Filter):
filter_type = 'particleproduction'

2
openmc/tallies.py
View file

@ -53,7 +53,7 @@ _NUCLIDE_CLASSES = (str, CrossNuclide, AggregateNuclide)
_FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
# Valid types of estimators
ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'}
ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog', 'next-event'}
class Tally(IDManagerMixin):

50
src/angle_energy.cpp Normal file
View file

@ -0,0 +1,50 @@
#include "openmc/angle_energy.h"
#include <algorithm> // for clamp
#include <cmath> // for sqrt
#include "openmc/constants.h"
#include "openmc/random_lcg.h"
namespace openmc {
double get_jac_and_transform(
double E_in, double& mu, double& E_out, uint64_t* seed, double awr)
{
double E_com = E_in / ((awr + 1.0) * (awr + 1.0));
return get_jac_and_transform_impl(E_com, mu, E_out, seed, awr);
}
double get_jac_and_transform_impl(
double E_com, double& mu, double& E_out, uint64_t* seed, double awr)
{
double E_cm = E_out;
double mu_lab = mu;
double D = mu_lab * mu_lab - 1.0 + E_cm / E_com;
if (D <= 0.0)
return 0.0;
D = std::sqrt(D);
if ((mu_lab <= 0.0) && (E_cm <= E_com))
return 0.0;
double E_out1 = E_com * (mu_lab + D) * (mu_lab + D);
double mult;
if (E_cm > E_com) {
mult = 1.0;
E_out = E_out1;
} else {
mult = 2.0;
if (prn(seed) < 0.5) {
E_out = E_out1;
} else {
E_out = E_com * (mu_lab - D) * (mu_lab - D);
}
}
mu = mu_lab * std::sqrt(E_out / E_cm) - std::sqrt(E_com / E_cm);
if ((std::abs(mu) > 1.0) && (std::abs(mu) < 1.0 + FP_PRECISION))
mu = std::clamp(mu, -1.0, 1.0);
return mult * E_out / (D * std::sqrt(E_cm * E_com));
}
} // namespace openmc

8
src/chain.cpp
View file

@ -74,9 +74,13 @@ void DecayPhotonAngleEnergy::sample(
mu = Uniform(-1., 1.).sample(seed).first;
}
double DecayPhotonAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double DecayPhotonAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
if (is_com)
fatal_error(
"DecayPhotonAngleEnergy distribution must be in lab coordinates");
E_out = photon_energy_->sample(seed).first;
return 0.5;
}

38
src/physics.cpp
View file

@ -24,7 +24,9 @@
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/next_event_scoring.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/thermal.h"
#include "openmc/weight_windows.h"
@ -156,6 +158,8 @@ void sample_neutron_reaction(Particle& p)
// exiting neutron
const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat();
if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) {
if (!model::active_point_tallies.empty())
fatal_error("Next-Event estimator does not support ncrystal materials");
ncrystal_mat.scatter(p);
} else {
scatter(p, i_nuclide);
@ -743,7 +747,7 @@ void scatter(Particle& p, int i_nuclide)
// Perform collision physics for inelastic scattering
const auto& rx {nuc->reactions_[i]};
inelastic_scatter(*nuc, *rx, p);
inelastic_scatter(i_nuclide, *rx, p);
p.event_mt() = rx->mt_;
}
@ -789,6 +793,10 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p)
// Find speed of neutron in CM
vel = v_n.norm();
if (!model::active_point_tallies.empty()) {
score_point_tally_elastic(p, i_nuclide, rx, 0, v_t);
}
// Sample scattering angle, checking if angle distribution is present (assume
// isotropic otherwise)
double mu_cm;
@ -836,8 +844,13 @@ void sab_scatter(int i_nuclide, int i_sab, Particle& p)
// Sample energy and angle
double E_out;
data::thermal_scatt[i_sab]->data_[i_temp].sample(
micro, p.E(), &E_out, &p.mu(), p.current_seed());
auto& sab = data::thermal_scatt[i_sab]->data_[i_temp];
if (!model::active_point_tallies.empty()) {
score_point_tally_sab(p, i_nuclide, sab, micro);
}
sab.sample(micro, p.E(), &E_out, &p.mu(), p.current_seed());
// Set energy to outgoing, change direction of particle
p.E() = E_out;
@ -1102,6 +1115,10 @@ void sample_fission_neutron(
site->delayed_group = 0;
}
if (!model::active_point_tallies.empty()) {
score_point_tally_fission(p, i_nuclide, rx, site->delayed_group);
}
// sample from prompt neutron energy distribution
int n_sample = 0;
double mu;
@ -1127,8 +1144,11 @@ void sample_fission_neutron(
site->u = rotate_angle(p.u(), mu, nullptr, seed);
}
void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p)
{
// Get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide]};
// copy energy of neutron
double E_in = p.E();
@ -1137,13 +1157,19 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
double mu;
rx.products_[0].sample(E_in, E, mu, p.current_seed());
double yield = (*rx.products_[0].yield_)(E_in);
if (!model::active_point_tallies.empty()) {
score_point_tally_inelastic(p, i_nuclide, rx, 0, yield);
}
// if scattering system is in center-of-mass, transfer cosine of scattering
// angle and outgoing energy from CM to LAB
if (rx.scatter_in_cm_) {
double E_cm = E;
// determine outgoing energy in lab
double A = nuc.awr_;
double A = nuc->awr_;
E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) /
((A + 1.0) * (A + 1.0));
@ -1164,8 +1190,6 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
// change direction of particle
p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
// evaluate yield
double yield = (*rx.products_[0].yield_)(E_in);
if (std::floor(yield) == yield && yield > 0) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {

65
src/ray.cpp
View file

@ -2,6 +2,8 @@
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
#include "openmc/mgxs_interface.h"
#include "openmc/settings.h"
namespace openmc {
@ -11,7 +13,7 @@ void Ray::compute_distance()
boundary() = distance_to_boundary(*this);
}
void Ray::trace()
void Ray::trace(double max_distance)
{
// To trace the ray from its origin all the way through the model, we have
// to proceed in two phases. In the first, the ray may or may not be found
@ -24,6 +26,8 @@ void Ray::trace()
// After phase one is done, we can starting tracing from cell to cell within
// the model. This step can use neighbor lists to accelerate the ray tracing.
double max = max_distance;
bool inside_cell;
// Check for location if the particle is already known
if (lowest_coord().cell() == C_NONE) {
@ -118,10 +122,20 @@ void Ray::trace()
// distance to properly check cell inclusion.
boundary().distance() += TINY_BIT;
double distance = std::min(boundary().distance(), max);
// Advance particle, prepare for next intersection
for (int lev = 0; lev < n_coord(); ++lev) {
coord(lev).r() += boundary().distance() * coord(lev).u();
coord(lev).r() += distance * coord(lev).u();
}
max -= distance;
if (max == 0.0) {
update_distance();
break;
}
surface() = boundary().surface();
// Initialize last cells from the current cell, because the cell() variable
// does not contain the data for the case of a single-segment ray
@ -136,8 +150,8 @@ void Ray::trace()
cross_lattice(*this, boundary(), settings::verbosity >= 10);
}
// Record how far the ray has traveled
traversal_distance_ += boundary().distance();
update_distance();
inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
// Call the specialized logic for this type of ray. Note that we do not
@ -165,4 +179,47 @@ void Ray::trace()
}
}
void Ray::update_distance()
{
// Record how far the ray has traveled
traversal_distance_ += boundary().distance();
}
void ParticleRay::on_intersection() {}
void ParticleRay::update_distance()
{
Ray::update_distance();
time() += boundary().distance() / speed();
// Calculate microscopic and macroscopic cross sections
if (material() != MATERIAL_VOID) {
if (settings::run_CE) {
if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
density_mult() != density_mult_last()) {
// If the material is the same as the last material and the
// temperature hasn't changed, we don't need to lookup cross
// sections again.
model::materials[material()]->calculate_xs(*this);
}
} else {
// Get the MG data; unlike the CE case above, we have to re-calculate
// cross sections for every collision since the cross sections may
// be angle-dependent
data::mg.macro_xs_[material()].calculate_xs(*this);
// Update the particle's group while we know we are multi-group
g_last() = g();
}
} else {
macro_xs().total = 0.0;
macro_xs().absorption = 0.0;
macro_xs().fission = 0.0;
macro_xs().nu_fission = 0.0;
}
traversal_mfp_ += macro_xs().total * boundary().distance();
}
} // namespace openmc

7
src/reaction_product.cpp
View file

@ -134,10 +134,11 @@ void ReactionProduct::sample(
sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
}
double ReactionProduct::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double ReactionProduct::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
return sample_dist(E_in, seed)
.sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
}
} // namespace openmc

11
src/secondary_correlated.cpp
View file

@ -260,10 +260,15 @@ void CorrelatedAngleEnergy::sample(
mu = sample_dist(E_in, E_out, seed).sample(seed).first;
}
double CorrelatedAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double CorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
return sample_dist(E_in, E_out, seed).evaluate(mu);
auto& dist = sample_dist(E_in, E_out, seed);
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
return jac * dist.evaluate(mu);
}
} // namespace openmc

10
src/secondary_kalbach.cpp
View file

@ -232,14 +232,18 @@ void KalbachMann::sample(
mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
}
}
double KalbachMann::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double KalbachMann::sample_energy_and_pdf(double E_in, double mu, double& E_out,
uint64_t* seed, bool is_com, double awr) const
{
double km_r, km_a;
sample_params(E_in, E_out, km_a, km_r, seed);
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
// https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
return km_a / (2 * std::sinh(km_a)) *
return jac * km_a / (2 * std::sinh(km_a)) *
(std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
}

9
src/secondary_nbody.cpp
View file

@ -72,11 +72,14 @@ void NBodyPhaseSpace::sample(
E_out = sample_energy(E_in, seed);
}
double NBodyPhaseSpace::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double NBodyPhaseSpace::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
E_out = sample_energy(E_in, seed);
return 0.5;
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
return jac * 0.5;
}
} // namespace openmc

68
src/secondary_thermal.cpp
View file

@ -53,11 +53,16 @@ void CoherentElasticAE::sample(
mu = 1.0 - 2.0 * energies[k] / E_in;
}
double CoherentElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double CoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
E_out = E_in;
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
const auto& factors = xs_.factors();
if (E_in < bragg_edges_.front())
@ -68,8 +73,9 @@ double CoherentElasticAE::sample_energy_and_pdf(
double E = 0.5 * (1 - mu) * E_in;
double C = 0.5 * E_in / factors[i];
return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
return jac * C *
get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
}
//==============================================================================
@ -90,15 +96,20 @@ void IncoherentElasticAE::sample(
double c = 2 * E_in * debye_waller_;
mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
}
double IncoherentElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double IncoherentElasticAE::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
E_out = E_in;
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
double c = 2 * E_in * debye_waller_;
double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
return A * std::exp(-c * (1 - mu));
return jac * A * std::exp(-c * (1 - mu));
}
//==============================================================================
@ -155,8 +166,8 @@ void IncoherentElasticAEDiscrete::sample(
E_out = E_in;
}
double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double IncoherentElasticAEDiscrete::sample_energy_and_pdf(double E_in,
double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
{
// Get index and interpolation factor for elastic grid
int i;
@ -165,8 +176,12 @@ double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
// Energy doesn't change in elastic scattering
E_out = E_in;
return get_pdf_discrete_interpolated(
mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu);
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, tensor::all),
mu_out_.slice(i + 1, tensor::all), f, mu);
}
//==============================================================================
@ -253,8 +268,8 @@ void IncoherentInelasticAEDiscrete::sample(
mu = (1 - f) * mu_ijk + f * mu_i1jk;
}
double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(double E_in,
double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const
{
// Get index and interpolation factor for inelastic grid
int i;
@ -263,8 +278,12 @@ double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
int j;
sample_params(E_in, E_out, j, seed);
return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
mu_out_.slice(i + 1, j, tensor::all), f, mu);
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
mu_out_.slice(i + 1, j, tensor::all), f, mu);
}
//==============================================================================
@ -403,8 +422,8 @@ void IncoherentInelasticAE::sample(
mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
}
double IncoherentInelasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double IncoherentInelasticAE::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
double f;
int l, j;
@ -412,8 +431,12 @@ double IncoherentInelasticAE::sample_energy_and_pdf(
const auto& mu_l = distribution_[l].mu;
return get_pdf_discrete_interpolated(
mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu);
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
return jac * get_pdf_discrete_interpolated(mu_l.slice(j, tensor::all),
mu_l.slice(j + 1, tensor::all), f, mu);
}
//==============================================================================
@ -458,10 +481,11 @@ void MixedElasticAE::sample(
sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
}
double MixedElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double MixedElasticAE::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
return sample_dist(E_in, seed)
.sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
}
} // namespace openmc

12
src/secondary_uncorrelated.cpp
View file

@ -65,8 +65,8 @@ void UncorrelatedAngleEnergy::sample(
E_out = energy_->sample(E_in, seed);
}
double UncorrelatedAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
double UncorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
double& E_out, uint64_t* seed, bool is_com, double awr) const
{
// Sample outgoing energy
if (energy_ != nullptr) {
@ -75,11 +75,15 @@ double UncorrelatedAngleEnergy::sample_energy_and_pdf(
E_out = E_in;
}
double jac = 1.0;
if (is_com)
jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);
if (!angle_.empty()) {
return angle_.evaluate(E_in, mu);
return jac * angle_.evaluate(E_in, mu);
} else {
// no angle distribution given => assume isotropic for all energies
return 0.5;
return jac * 0.5;
}
}

1
src/simulation.cpp
View file

@ -339,6 +339,7 @@ double k_abs_tra {0.0};
double log_spacing;
int n_lost_particles {0};
bool need_depletion_rx {false};
bool nonvacuum_boundary_present {false};
int restart_batch;
bool satisfy_triggers {false};
int ssw_current_file;

6
src/source.cpp
View file

@ -38,6 +38,8 @@
#include "openmc/simulation.h"
#include "openmc/state_point.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/next_event_scoring.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -1207,6 +1209,10 @@ SourceSite sample_external_source(uint64_t* seed)
site.E = data::mg.num_energy_groups_ - site.E - 1.;
}
if (!model::active_point_tallies.empty()) {
score_point_tally_source(site, i);
}
return site;
}

2
src/state_point.cpp
View file

@ -221,6 +221,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
write_dataset(tally_group, "estimator", "tracklength");
} else if (tally->estimator_ == TallyEstimator::COLLISION) {
write_dataset(tally_group, "estimator", "collision");
} else if (tally->estimator_ == TallyEstimator::NEXT_EVENT) {
write_dataset(tally_group, "estimator", "next-event");
}
write_dataset(tally_group, "n_realizations", tally->n_realizations_);

4
src/surface.cpp
View file

@ -17,6 +17,7 @@
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
@ -1175,6 +1176,7 @@ void read_surfaces(pugi::xml_node node,
std::unordered_map<int, double>& albedo_map,
std::unordered_map<int, int>& periodic_sense_map)
{
simulation::nonvacuum_boundary_present = false;
// Count the number of surfaces
int n_surfaces = 0;
for (pugi::xml_node surf_node : node.children("surface")) {
@ -1245,6 +1247,8 @@ void read_surfaces(pugi::xml_node node,
// Check for a periodic surface
if (check_for_node(surf_node, "boundary")) {
std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
if (surf_bc != "vacuum")
simulation::nonvacuum_boundary_present = true;
if (surf_bc == "periodic") {
periodic_sense_map[model::surfaces.back()->id_] = 0;
// Check for surface albedo. Skip sanity check as it is already done

3
src/tallies/filter.cpp
View file

@ -32,6 +32,7 @@
#include "openmc/tallies/filter_parent_nuclide.h"
#include "openmc/tallies/filter_particle.h"
#include "openmc/tallies/filter_particle_production.h"
#include "openmc/tallies/filter_point.h"
#include "openmc/tallies/filter_polar.h"
#include "openmc/tallies/filter_reaction.h"
#include "openmc/tallies/filter_sph_harm.h"
@ -150,6 +151,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
return Filter::create<ParticleFilter>(id);
} else if (type == "particleproduction") {
return Filter::create<ParticleProductionFilter>(id);
} else if (type == "point") {
return Filter::create<PointFilter>(id);
} else if (type == "polar") {
return Filter::create<PolarFilter>(id);
} else if (type == "reaction") {

79
src/tallies/filter_point.cpp Normal file
View file

@ -0,0 +1,79 @@
#include "openmc/tallies/filter_point.h"
#include "openmc/tallies/tally_scoring.h"
#include <fmt/core.h>
#include "openmc/math_functions.h"
#include "openmc/simulation.h"
#include "openmc/xml_interface.h"
namespace openmc {
void PointFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(node, "bins");
// Convert to vector of detectors
vector<std::pair<Position, double>> detectors;
size_t n = bins.size() / 4;
for (int i = 0; i < n; ++i) {
Position pos {bins[4 * i], bins[4 * i + 1], bins[4 * i + 2]};
detectors.push_back(std::make_pair(pos, bins[4 * i + 3]));
}
this->set_detectors(detectors);
}
void PointFilter::set_detectors(span<std::pair<Position, double>> detectors)
{
// Clear existing detectors
detectors_.clear();
detectors_.reserve(detectors.size());
// Set detectors and number of bins
for (auto d : detectors) {
detectors_.push_back(d);
}
n_bins_ = detectors_.size();
}
void PointFilter::get_all_bins(
const Particle& p, TallyEstimator estimator, FilterMatch& match) const
{
double attenuation = p.wgt() / p.wgt_last();
int i = 0;
for (auto [pos, r] : detectors_) {
if ((p.r() - pos).norm() < FP_COINCIDENT) {
match.bins_.push_back(i);
double weight;
double distance = (p.r_last() - pos).norm();
if (distance > r) {
weight = attenuation / (distance * distance);
} else {
weight = 3.0 * exprel(-p.macro_xs().total * r) / (r * r);
}
match.weights_.push_back(weight);
}
++i;
}
}
void PointFilter::to_statepoint(hid_t filter_group) const
{
Filter::to_statepoint(filter_group);
vector<double> detectors;
for (auto [pos, r] : detectors_) {
detectors.push_back(pos[0]);
detectors.push_back(pos[1]);
detectors.push_back(pos[2]);
detectors.push_back(r);
}
write_dataset(filter_group, "bins", detectors);
}
std::string PointFilter::text_label(int bin) const
{
auto [pos, r] = detectors_.at(bin);
return fmt::format("Point: {} {}", pos, r);
}
} // namespace openmc

115
src/tallies/next_event_scoring.cpp Normal file
View file

@ -0,0 +1,115 @@
#include "openmc/tallies/next_event_scoring.h"
#include <algorithm> // for clamp
#include "openmc/secondary_uncorrelated.h"
#include "openmc/source.h"
namespace openmc {
void score_point_tally_elastic(
Particle& p, int i_nuclide, const Reaction& rx, int i_product, Direction v_t)
{
const auto& nuc {data::nuclides[i_nuclide]};
double awr = nuc->awr_;
// Neutron velocity in LAB
Direction v_n = std::sqrt(p.E()) * p.u();
auto u_n = v_n / v_n.norm();
// Velocity of center-of-mass
Direction v_cm = (v_n + awr * v_t) / (awr + 1.0);
auto u_cm = v_cm / v_cm.norm();
double E_in = p.E();
double E_com = v_cm.dot(v_cm);
double E_out = (v_n - v_cm).dot(v_n - v_cm);
auto& d = rx.products_[i_product].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
auto pdf = [&](Direction u, double& E) {
double mu = u.dot(u_cm);
double mu_l = u.dot(u_n);
E = E_out;
double jac =
get_jac_and_transform_impl(E_com, mu, E, p.current_seed(), awr);
double mu_cm =
1.0 + mu_l * std::sqrt(E_in * E) / E_out - (E_in + E) / (2.0 * E_out);
mu_cm = std::clamp(mu_cm, -1.0, 1.0);
if (!d_->angle().empty()) {
return jac * d_->angle().evaluate(p.E(), mu_cm) / (2.0 * PI);
} else {
return jac * 0.5 / (2.0 * PI);
}
};
score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
}
void score_point_tally_inelastic(
Particle& p, int i_nuclide, const Reaction& rx, int i_product, double yield)
{
const auto& nuc {data::nuclides[i_nuclide]};
double awr = nuc->awr_;
auto u_n = p.u();
auto E_n = p.E();
auto is_com = rx.scatter_in_cm_;
auto pdf = [&](Direction u, double& E) {
double mu = u.dot(u_n);
return rx.products_[i_product].sample_energy_and_pdf(
E_n, mu, E, p.current_seed(), is_com, awr) /
(2.0 * PI) * yield;
};
score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
}
void score_point_tally_fission(
Particle& p, int i_nuclide, const Reaction& rx, int i_product)
{
const auto& nuc {data::nuclides[i_nuclide]};
double awr = nuc->awr_;
auto u_n = p.u();
auto E_n = p.E();
auto is_com = rx.scatter_in_cm_;
auto pdf = [&](Direction u, double& E) {
double mu = u.dot(u_n);
return rx.products_[i_product].sample_energy_and_pdf(
E_n, mu, E, p.current_seed(), is_com, awr) /
(2.0 * PI);
};
score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
}
void score_point_tally_sab(Particle& p, int i_nuclide, const ThermalData& sab,
const NuclideMicroXS& micro)
{
const auto& nuc {data::nuclides[i_nuclide]};
double awr = nuc->awr_;
auto u_n = p.u();
auto E_n = p.E();
auto pdf = [&](Direction u, double& E) {
double mu = u.dot(u_n);
return sab.sample_energy_and_pdf(
micro, E_n, mu, E, p.current_seed(), false, awr) /
(2.0 * PI);
};
score_point_tally_impl(p.r(), p.type(), p.time(), pdf);
}
void score_point_tally_source(SourceSite& site, int source_index)
{
auto src_ = model::external_sources[source_index].get();
auto src = dynamic_cast<IndependentSource*>(src_);
if (!src)
fatal_error("Only independent source is valid for point detectors.");
auto pdf = [&](Direction u, double& E) {
E = site.E;
return src->angle()->evaluate(u);
};
score_point_tally_impl(site.r, site.particle, site.time, pdf);
}
} // namespace openmc

48
src/tallies/tally.cpp
View file

@ -31,6 +31,7 @@
#include "openmc/tallies/filter_meshmaterial.h"
#include "openmc/tallies/filter_meshsurface.h"
#include "openmc/tallies/filter_particle.h"
#include "openmc/tallies/filter_point.h"
#include "openmc/tallies/filter_sph_harm.h"
#include "openmc/tallies/filter_surface.h"
#include "openmc/tallies/filter_time.h"
@ -60,11 +61,13 @@ vector<int> active_analog_tallies;
vector<int> active_tracklength_tallies;
vector<int> active_timed_tracklength_tallies;
vector<int> active_collision_tallies;
vector<int> active_point_tallies;
vector<int> active_meshsurf_tallies;
vector<int> active_surface_tallies;
vector<int> active_pulse_height_tallies;
vector<int32_t> pulse_height_cells;
vector<double> time_grid;
std::set<Position> active_point_detectors;
} // namespace model
namespace simulation {
@ -540,6 +543,7 @@ void Tally::set_scores(const vector<std::string>& scores)
bool legendre_present = false;
bool cell_present = false;
bool cellfrom_present = false;
bool point_present = false;
bool material_present = false;
bool materialfrom_present = false;
bool surface_present = false;
@ -562,6 +566,10 @@ void Tally::set_scores(const vector<std::string>& scores)
materialfrom_present = true;
} else if (filt->type() == FilterType::MATERIAL) {
material_present = true;
} else if (filt->type() == FilterType::POINT) {
point_present = true;
type_ = TallyType::POINT;
estimator_ = TallyEstimator::NEXT_EVENT;
} else if (filt->type() == FilterType::SURFACE) {
surface_present = true;
} else if (filt->type() == FilterType::MESH_SURFACE) {
@ -574,6 +582,21 @@ void Tally::set_scores(const vector<std::string>& scores)
(surface_present || cell_present || cellfrom_present || material_present ||
materialfrom_present);
if (point_present) {
if (!settings::run_CE)
fatal_error("Cannot use point detectors in multi-group mode.");
if (simulation::nonvacuum_boundary_present)
fatal_error(
"Cannot use point detectors with non-vacuum boundary conditions.");
if (legendre_present)
fatal_error("Cannot use LegendreFilter with PointFilter.");
if (energyout_present)
fatal_error("Cannot use EnergyoutFilter with PointFilter.");
if (surface_present || meshsurface_present)
fatal_error(
"Cannot use surface or mesh-surface filters with PointFilter.");
}
// Iterate over the given scores.
for (auto score_str : scores) {
// Make sure a delayed group filter wasn't used with an incompatible
@ -623,6 +646,9 @@ void Tally::set_scores(const vector<std::string>& scores)
case SCORE_NU_SCATTER:
if (settings::run_CE) {
if (point_present)
fatal_error("Cannot use nu-scatter score with PointFilter in "
"continuous energy mode.");
estimator_ = TallyEstimator::ANALOG;
} else {
if (energyout_present || legendre_present)
@ -631,6 +657,8 @@ void Tally::set_scores(const vector<std::string>& scores)
break;
case SCORE_CURRENT:
if (point_present)
fatal_error("Cannot use current score with PointFilter.");
// Check which type of current is desired: mesh or surface currents.
if (meshsurface_present) {
if (non_meshsurface_types_present)
@ -644,6 +672,8 @@ void Tally::set_scores(const vector<std::string>& scores)
break;
case HEATING:
if (point_present)
fatal_error("Cannot use heating score with PointFilter.");
if (settings::photon_transport) {
// Photon heating requires a collision estimator (analog energy
// balance). However, if the tally only scores neutrons, we can keep the
@ -665,6 +695,8 @@ void Tally::set_scores(const vector<std::string>& scores)
break;
case SCORE_PULSE_HEIGHT: {
if (point_present)
fatal_error("Cannot use pulse-height score with PointFilter.");
if (non_cell_energy_present) {
fatal_error("Pulse-height tallies are not compatible with filters "
"other than CellFilter and EnergyFilter");
@ -686,6 +718,8 @@ void Tally::set_scores(const vector<std::string>& scores)
case SCORE_IFP_TIME_NUM:
case SCORE_IFP_BETA_NUM:
case SCORE_IFP_DENOM:
if (point_present)
fatal_error("Cannot use ifp scores with PointFilter.");
estimator_ = TallyEstimator::COLLISION;
break;
}
@ -1184,6 +1218,8 @@ void setup_active_tallies()
model::active_meshsurf_tallies.clear();
model::active_surface_tallies.clear();
model::active_pulse_height_tallies.clear();
model::active_point_tallies.clear();
model::active_point_detectors.clear();
model::time_grid.clear();
for (auto i = 0; i < model::tallies.size(); ++i) {
@ -1225,6 +1261,16 @@ void setup_active_tallies()
case TallyType::PULSE_HEIGHT:
model::active_pulse_height_tallies.push_back(i);
break;
case TallyType::POINT:
model::active_point_tallies.push_back(i);
// Populate the set of unique detector positions from PointFilter
if (auto pf = tally.get_filter<PointFilter>()) {
for (const auto& [pos, r0] : pf->detectors()) {
model::active_point_detectors.insert(pos);
}
}
break;
}
}
}
@ -1248,6 +1294,8 @@ void free_memory_tally()
model::active_meshsurf_tallies.clear();
model::active_surface_tallies.clear();
model::active_pulse_height_tallies.clear();
model::active_point_tallies.clear();
model::active_point_detectors.clear();
model::time_grid.clear();
model::tally_map.clear();

3
src/tallies/tally_scoring.cpp
View file

@ -11,8 +11,10 @@
#include "openmc/photon.h"
#include "openmc/reaction_product.h"
#include "openmc/search.h"
#include "openmc/secondary_uncorrelated.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/string_utils.h"
#include "openmc/surface.h"
#include "openmc/tallies/derivative.h"
@ -2785,4 +2787,5 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
p.coord(0).cell() = orig_cell;
p.E_last() = orig_E_last;
}
} // namespace openmc

5
src/thermal.cpp
View file

@ -314,10 +314,11 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
}
double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
double E_in, double mu, double& E_out, uint64_t* seed) const
double E_in, double mu, double& E_out, uint64_t* seed, bool is_com,
double awr) const
{
return sample_dist(micro_xs, E_in, seed)
.sample_energy_and_pdf(E_in, mu, E_out, seed);
.sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
}
void free_memory_thermal()