17271733860

Committed 27 Aug 2025 03:46PM UTC coverage: 84.486% (-0.7%) from 85.204%

Build # 17271733860

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge 06ae298be into d1df80a21

Pull Request Pull Request #3550: [Point Detector] Add distribution get_pdf functionality

Run Details

4 of 551 new or added lines in 11 files covered. (0.73%)

4 existing lines in 2 files now uncovered.

52940 of 62661 relevant lines covered (84.49%)

37443356.37 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

65.56

/src/reaction_product.cpp

#include "openmc/reaction_product.h"

#include <string> // for string

#include <fmt/core.h>

#include "openmc/endf.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/memory.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/secondary_correlated.h"
#include "openmc/secondary_kalbach.h"
#include "openmc/secondary_nbody.h"
#include "openmc/secondary_thermal.h"
#include "openmc/secondary_uncorrelated.h"
#include "openmc/tallies/tally_scoring.h"

namespace openmc {

//==============================================================================
// ReactionProduct implementation
//==============================================================================

ReactionProduct::ReactionProduct(hid_t group)
{
  // Read particle type
  std::string temp;
  read_attribute(group, "particle", temp);
  particle_ = str_to_particle_type(temp);

  // Read emission mode and decay rate
  read_attribute(group, "emission_mode", temp);
  if (temp == "prompt") {
    emission_mode_ = EmissionMode::prompt;
  } else if (temp == "delayed") {
    emission_mode_ = EmissionMode::delayed;
  } else if (temp == "total") {
    emission_mode_ = EmissionMode::total;
  }

  // Read decay rate for delayed emission
  if (emission_mode_ == EmissionMode::delayed) {
    if (attribute_exists(group, "decay_rate")) {
      read_attribute(group, "decay_rate", decay_rate_);
    } else if (particle_ == ParticleType::neutron) {
      warning(fmt::format("Decay rate doesn't exist for delayed neutron "
                          "emission ({}).",
        object_name(group)));
    }
  }

  // Read secondary particle yield
  yield_ = read_function(group, "yield");

  int n;
  read_attribute(group, "n_distribution", n);

  for (int i = 0; i < n; ++i) {
    std::string s {"distribution_"};
    s.append(std::to_string(i));
    hid_t dgroup = open_group(group, s.c_str());

    // Read applicability
    if (n > 1) {
      hid_t app = open_dataset(dgroup, "applicability");
      applicability_.emplace_back(app);
      close_dataset(app);
    }

    // Determine distribution type and read data
    read_attribute(dgroup, "type", temp);
    if (temp == "uncorrelated") {
      distribution_.push_back(make_unique<UncorrelatedAngleEnergy>(dgroup));
    } else if (temp == "correlated") {
      distribution_.push_back(make_unique<CorrelatedAngleEnergy>(dgroup));
    } else if (temp == "nbody") {
      distribution_.push_back(make_unique<NBodyPhaseSpace>(dgroup));
    } else if (temp == "kalbach-mann") {
      distribution_.push_back(make_unique<KalbachMann>(dgroup));
    }

    close_group(dgroup);
  }
}

ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
{
  particle_ = ParticleType::photon;
  emission_mode_ = EmissionMode::delayed;

  // Get chain nuclide object for radionuclide
  parent_nuclide_ = data::chain_nuclide_map.at(product.name);
  const auto& chain_nuc = data::chain_nuclides[parent_nuclide_].get();

  // Determine decay constant in [s^-1]
  decay_rate_ = chain_nuc->decay_constant();

  // Determine number of photons per decay and set yield
  double photon_per_sec = chain_nuc->photon_energy()->integral();
  double photon_per_decay = photon_per_sec / decay_rate_;
  vector<double> coef = {product.branching_ratio * photon_per_decay};
  yield_ = make_unique<Polynomial>(coef);

  // Set decay photon angle-energy distribution
  distribution_.push_back(
    make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));
}

void ReactionProduct::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  auto n = applicability_.size();
  if (n > 1) {
    double prob = 0.0;
    double c = prn(seed);
    for (int i = 0; i < n; ++i) {
      // Determine probability that i-th energy distribution is sampled
      prob += applicability_[i](E_in);

      // If i-th distribution is sampled, sample energy from the distribution
      if (c <= prob) {
        distribution_[i]->sample(E_in, E_out, mu, seed);
        break;
      }
    }
  } else {
    // If only one distribution is present, go ahead and sample it
    distribution_[0]->sample(E_in, E_out, mu, seed);
  }
}
void ReactionProduct::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
{
  /* function to get detector position from user - Future implementation
   double det_pos[4];
   get_det_pos(det_pos, i_tally);
  */
  double det_pos[4] = {0.0, 0.0, 0.0, 1.0}; // Placeholder for detector position

  int distribution_index;
  auto n = applicability_.size();
  if (n > 1) {
    double prob = 0.0;
    double c = prn(seed);
    for (int i = 0; i < n; ++i) {
      // Determine probability that i-th energy distribution is sampled
      prob += applicability_[i](E_in);

      // If i-th distribution is sampled, sample energy from the distribution
      if (c <= prob) {
        // distribution_[i]->sample(E_in, E_out, mu, seed);
        distribution_index = i;
        break;
      }
    }
  } else {
    // If only one distribution is present, go ahead and sample it
    // distribution_[0]->sample(E_in, E_out, mu, seed);
    distribution_index = 0;
  }
  // now extract pdf

  AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();

  if (CorrelatedAngleEnergy* correlatedAE =
        dynamic_cast<CorrelatedAngleEnergy*>(angleEnergyPtr)) {

    (*correlatedAE)
      .get_pdf(
        det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
    // Handle CorrelatedAngleEnergy
  } else if (KalbachMann* kalbachMann =
               dynamic_cast<KalbachMann*>(angleEnergyPtr)) {

    (*kalbachMann)
      .get_pdf(
        det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);

    // Handle KalbachMann
  } else if (NBodyPhaseSpace* nBodyPS =
               dynamic_cast<NBodyPhaseSpace*>(angleEnergyPtr)) {

    (*nBodyPS).get_pdf(
      det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
    // Handle NBodyPhaseSpace
  } else if (UncorrelatedAngleEnergy* uncorrelatedAE =
               dynamic_cast<UncorrelatedAngleEnergy*>(angleEnergyPtr)) {

    (*uncorrelatedAE)
      .get_pdf(
        det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
    // Handle UncorrelatedAngleEnergy
  } else {
    std::cout << "Unknown derived type." << std::endl;
  }
}

} // namespace openmc

1	#include "openmc/reaction_product.h"
2
3	#include <string> // for string
4
5	#include <fmt/core.h>
6
7	#include "openmc/endf.h"
8	#include "openmc/error.h"
9	#include "openmc/hdf5_interface.h"
10	#include "openmc/memory.h"
11	#include "openmc/particle.h"
12	#include "openmc/random_lcg.h"
13	#include "openmc/secondary_correlated.h"
14	#include "openmc/secondary_kalbach.h"
15	#include "openmc/secondary_nbody.h"
16	#include "openmc/secondary_thermal.h"
17	#include "openmc/secondary_uncorrelated.h"
18	#include "openmc/tallies/tally_scoring.h"
19
20	namespace openmc {
21
22	//==============================================================================
23	// ReactionProduct implementation
24	//==============================================================================
25
26	ReactionProduct::ReactionProduct(hid_t group)	3,008,627✔
27	{
28	// Read particle type
29	std::string temp;	3,008,627✔
30	read_attribute(group, "particle", temp);	3,008,627✔
31	particle_ = str_to_particle_type(temp);	3,008,627✔
32
33	// Read emission mode and decay rate
34	read_attribute(group, "emission_mode", temp);	3,008,627✔
35	if (temp == "prompt") {	3,008,627✔
36	emission_mode_ = EmissionMode::prompt;	2,940,743✔
37	} else if (temp == "delayed") {	67,884✔
38	emission_mode_ = EmissionMode::delayed;	67,884✔
39	} else if (temp == "total") {	×
40	emission_mode_ = EmissionMode::total;	×
41	}
42
43	// Read decay rate for delayed emission
44	if (emission_mode_ == EmissionMode::delayed) {	3,008,627✔
45	if (attribute_exists(group, "decay_rate")) {	67,884✔
46	read_attribute(group, "decay_rate", decay_rate_);	67,884✔
47	} else if (particle_ == ParticleType::neutron) {	×
48	warning(fmt::format("Decay rate doesn't exist for delayed neutron "	×
49	"emission ({}).",
50	object_name(group)));	×
51	}
52	}
53
54	// Read secondary particle yield
55	yield_ = read_function(group, "yield");	3,008,627✔
56
57	int n;
58	read_attribute(group, "n_distribution", n);	3,008,627✔
59
60	for (int i = 0; i < n; ++i) {	6,018,179✔
61	std::string s {"distribution_"};	3,009,552✔
62	s.append(std::to_string(i));	3,009,552✔
63	hid_t dgroup = open_group(group, s.c_str());	3,009,552✔
64
65	// Read applicability
66	if (n > 1) {	3,009,552✔
67	hid_t app = open_dataset(dgroup, "applicability");	1,754✔
68	applicability_.emplace_back(app);	1,754✔
69	close_dataset(app);	1,754✔
70	}
71
72	// Determine distribution type and read data
73	read_attribute(dgroup, "type", temp);	3,009,552✔
74	if (temp == "uncorrelated") {	3,009,552✔
75	distribution_.push_back(make_unique<UncorrelatedAngleEnergy>(dgroup));	2,904,832✔
76	} else if (temp == "correlated") {	104,720✔
77	distribution_.push_back(make_unique<CorrelatedAngleEnergy>(dgroup));	44,030✔
78	} else if (temp == "nbody") {	60,690✔
79	distribution_.push_back(make_unique<NBodyPhaseSpace>(dgroup));	697✔
80	} else if (temp == "kalbach-mann") {	59,993✔
81	distribution_.push_back(make_unique<KalbachMann>(dgroup));	59,993✔
82	}
83
84	close_group(dgroup);	3,009,552✔
85	}	3,009,552✔
86	}	3,008,627✔
87
88	ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)	132✔
89	{
90	particle_ = ParticleType::photon;	132✔
91	emission_mode_ = EmissionMode::delayed;	132✔
92
93	// Get chain nuclide object for radionuclide
94	parent_nuclide_ = data::chain_nuclide_map.at(product.name);	132✔
95	const auto& chain_nuc = data::chain_nuclides[parent_nuclide_].get();	132✔
96
97	// Determine decay constant in [s^-1]
98	decay_rate_ = chain_nuc->decay_constant();	132✔
99
100	// Determine number of photons per decay and set yield
101	double photon_per_sec = chain_nuc->photon_energy()->integral();	132✔
102	double photon_per_decay = photon_per_sec / decay_rate_;	132✔
103	vector<double> coef = {product.branching_ratio * photon_per_decay};	132✔
104	yield_ = make_unique<Polynomial>(coef);	132✔
105
106	// Set decay photon angle-energy distribution
107	distribution_.push_back(	132✔
108	make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));	264✔
109	}	132✔
110
111	void ReactionProduct::sample(	36,552,002✔
112	double E_in, double& E_out, double& mu, uint64_t* seed) const
113	{
114	auto n = applicability_.size();	36,552,002✔
115	if (n > 1) {	36,552,002✔
116	double prob = 0.0;	1,078✔
117	double c = prn(seed);	1,078✔
118	for (int i = 0; i < n; ++i) {	1,881✔
119	// Determine probability that i-th energy distribution is sampled
120	prob += applicability_[i](E_in);	1,881✔
121
122	// If i-th distribution is sampled, sample energy from the distribution
123	if (c <= prob) {	1,881✔
124	distribution_[i]->sample(E_in, E_out, mu, seed);	1,078✔
125	break;	1,078✔
126	}
127	}
128	} else {
129	// If only one distribution is present, go ahead and sample it
130	distribution_[0]->sample(E_in, E_out, mu, seed);	36,550,924✔
131	}
132	}	36,552,002✔
NEW 133	void ReactionProduct::get_pdf(int i_tally, double E_in, double& E_out,	×
134	uint64_t* seed, Particle& p, std::vector<double>& mu_cm,
135	std::vector<double>& Js, std::vector<Particle>& ghost_particles,
136	std::vector<double>& pdfs_lab) const
137	{
138	/* function to get detector position from user - Future implementation
139	double det_pos[4];
140	get_det_pos(det_pos, i_tally);
141	*/
NEW 142	double det_pos[4] = {0.0, 0.0, 0.0, 1.0}; // Placeholder for detector position	×
143
144	int distribution_index;
NEW 145	auto n = applicability_.size();	×
NEW 146	if (n > 1) {	×
NEW 147	double prob = 0.0;	×
NEW 148	double c = prn(seed);	×
NEW 149	for (int i = 0; i < n; ++i) {	×
150	// Determine probability that i-th energy distribution is sampled
NEW 151	prob += applicability_[i](E_in);	×
152
153	// If i-th distribution is sampled, sample energy from the distribution
NEW 154	if (c <= prob) {	×
155	// distribution_[i]->sample(E_in, E_out, mu, seed);
NEW 156	distribution_index = i;	×
NEW 157	break;	×
158	}
159	}
160	} else {
161	// If only one distribution is present, go ahead and sample it
162	// distribution_[0]->sample(E_in, E_out, mu, seed);
NEW 163	distribution_index = 0;	×
164	}
165	// now extract pdf
166
NEW 167	AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();	×
168
NEW 169	if (CorrelatedAngleEnergy* correlatedAE =	×
NEW 170	dynamic_cast<CorrelatedAngleEnergy*>(angleEnergyPtr)) {	×
171
172	(*correlatedAE)
NEW 173	.get_pdf(	×
174	det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
175	// Handle CorrelatedAngleEnergy
NEW 176	} else if (KalbachMann* kalbachMann =	×
NEW 177	dynamic_cast<KalbachMann*>(angleEnergyPtr)) {	×
178
179	(*kalbachMann)
NEW 180	.get_pdf(	×
181	det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
182
183	// Handle KalbachMann
NEW 184	} else if (NBodyPhaseSpace* nBodyPS =	×
NEW 185	dynamic_cast<NBodyPhaseSpace*>(angleEnergyPtr)) {	×
186
NEW 187	(*nBodyPS).get_pdf(	×
188	det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
189	// Handle NBodyPhaseSpace
NEW 190	} else if (UncorrelatedAngleEnergy* uncorrelatedAE =	×
NEW 191	dynamic_cast<UncorrelatedAngleEnergy*>(angleEnergyPtr)) {	×
192
193	(*uncorrelatedAE)
NEW 194	.get_pdf(	×
195	det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab);
196	// Handle UncorrelatedAngleEnergy
197	} else {
NEW 198	std::cout << "Unknown derived type." << std::endl;	×
199	}
200	}
201
202	} // namespace openmc

openmc-dev / openmc / 17271733860

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous