17475377863

Committed 04 Sep 2025 08:05PM UTC coverage: 84.707% (-0.5%) from 85.209%

Build # 17475377863

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge ff83fc4b9 into 591856472

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

Run Details

0 of 370 new or added lines in 10 files covered. (0.0%)

3 existing lines in 1 file now uncovered.

52947 of 62506 relevant lines covered (84.71%)

38165206.92 hits per line

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

76.62

/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_uncorrelated.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);
  }
}

double ReactionProduct::sample_energy_and_pdf(
  double E_in, double mu, double& E_out, uint64_t* seed) const
{

  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 {
    distribution_index = 0;
  }
  // now extract pdf

  AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();
  return angleEnergyPtr->sample_energy_and_pdf(E_in, mu, E_out, seed);
}

} // 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_uncorrelated.h"
17
18	namespace openmc {
19
20	//==============================================================================
21	// ReactionProduct implementation
22	//==============================================================================
23
24	ReactionProduct::ReactionProduct(hid_t group)	3,299,479✔
25	{
26	// Read particle type
27	std::string temp;	3,299,479✔
28	read_attribute(group, "particle", temp);	3,299,479✔
29	particle_ = str_to_particle_type(temp);	3,299,479✔
30
31	// Read emission mode and decay rate
32	read_attribute(group, "emission_mode", temp);	3,299,479✔
33	if (temp == "prompt") {	3,299,479✔
34	emission_mode_ = EmissionMode::prompt;	3,224,215✔
35	} else if (temp == "delayed") {	75,264✔
36	emission_mode_ = EmissionMode::delayed;	75,264✔
37	} else if (temp == "total") {	×
38	emission_mode_ = EmissionMode::total;	×
39	}
40
41	// Read decay rate for delayed emission
42	if (emission_mode_ == EmissionMode::delayed) {	3,299,479✔
43	if (attribute_exists(group, "decay_rate")) {	75,264✔
44	read_attribute(group, "decay_rate", decay_rate_);	75,264✔
45	} else if (particle_ == ParticleType::neutron) {	×
46	warning(fmt::format("Decay rate doesn't exist for delayed neutron "	×
47	"emission ({}).",
48	object_name(group)));	×
49	}
50	}
51
52	// Read secondary particle yield
53	yield_ = read_function(group, "yield");	3,299,479✔
54
55	int n;
56	read_attribute(group, "n_distribution", n);	3,299,479✔
57
58	for (int i = 0; i < n; ++i) {	6,599,957✔
59	std::string s {"distribution_"};	3,300,478✔
60	s.append(std::to_string(i));	3,300,478✔
61	hid_t dgroup = open_group(group, s.c_str());	3,300,478✔
62
63	// Read applicability
64	if (n > 1) {	3,300,478✔
65	hid_t app = open_dataset(dgroup, "applicability");	1,902✔
66	applicability_.emplace_back(app);	1,902✔
67	close_dataset(app);	1,902✔
68	}
69
70	// Determine distribution type and read data
71	read_attribute(dgroup, "type", temp);	3,300,478✔
72	if (temp == "uncorrelated") {	3,300,478✔
73	distribution_.push_back(make_unique<UncorrelatedAngleEnergy>(dgroup));	3,186,783✔
74	} else if (temp == "correlated") {	113,695✔
75	distribution_.push_back(make_unique<CorrelatedAngleEnergy>(dgroup));	47,799✔
76	} else if (temp == "nbody") {	65,896✔
77	distribution_.push_back(make_unique<NBodyPhaseSpace>(dgroup));	730✔
78	} else if (temp == "kalbach-mann") {	65,166✔
79	distribution_.push_back(make_unique<KalbachMann>(dgroup));	65,166✔
80	}
81
82	close_group(dgroup);	3,300,478✔
83	}	3,300,478✔
84	}	3,299,479✔
85
86	ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)	132✔
87	{
88	particle_ = ParticleType::photon;	132✔
89	emission_mode_ = EmissionMode::delayed;	132✔
90
91	// Get chain nuclide object for radionuclide
92	parent_nuclide_ = data::chain_nuclide_map.at(product.name);	132✔
93	const auto& chain_nuc = data::chain_nuclides[parent_nuclide_].get();	132✔
94
95	// Determine decay constant in [s^-1]
96	decay_rate_ = chain_nuc->decay_constant();	132✔
97
98	// Determine number of photons per decay and set yield
99	double photon_per_sec = chain_nuc->photon_energy()->integral();	132✔
100	double photon_per_decay = photon_per_sec / decay_rate_;	132✔
101	vector<double> coef = {product.branching_ratio * photon_per_decay};	132✔
102	yield_ = make_unique<Polynomial>(coef);	132✔
103
104	// Set decay photon angle-energy distribution
105	distribution_.push_back(	132✔
106	make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));	264✔
107	}	132✔
108
109	void ReactionProduct::sample(	39,684,582✔
110	double E_in, double& E_out, double& mu, uint64_t* seed) const
111	{
112	auto n = applicability_.size();	39,684,582✔
113	if (n > 1) {	39,684,582✔
114	double prob = 0.0;	1,078✔
115	double c = prn(seed);	1,078✔
116	for (int i = 0; i < n; ++i) {	1,881✔
117	// Determine probability that i-th energy distribution is sampled
118	prob += applicability_[i](E_in);	1,881✔
119
120	// If i-th distribution is sampled, sample energy from the distribution
121	if (c <= prob) {	1,881✔
122	distribution_[i]->sample(E_in, E_out, mu, seed);	1,078✔
123	break;	1,078✔
124	}
125	}
126	} else {
127	// If only one distribution is present, go ahead and sample it
128	distribution_[0]->sample(E_in, E_out, mu, seed);	39,683,504✔
129	}
130	}	39,684,582✔
131
NEW 132	double ReactionProduct::sample_energy_and_pdf(	×
133	double E_in, double mu, double& E_out, uint64_t* seed) const
134	{
135
136	int distribution_index;
NEW 137	auto n = applicability_.size();	×
NEW 138	if (n > 1) {	×
NEW 139	double prob = 0.0;	×
NEW 140	double c = prn(seed);	×
NEW 141	for (int i = 0; i < n; ++i) {	×
142	// Determine probability that i-th energy distribution is sampled
NEW 143	prob += applicability_[i](E_in);	×
144
145	// If i-th distribution is sampled, sample energy from the distribution
NEW 146	if (c <= prob) {	×
147	// distribution_[i]->sample(E_in, E_out, mu, seed);
NEW 148	distribution_index = i;	×
NEW 149	break;	×
150	}
151	}
152	} else {
NEW 153	distribution_index = 0;	×
154	}
155	// now extract pdf
156
NEW 157	AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get();	×
NEW 158	return angleEnergyPtr->sample_energy_and_pdf(E_in, mu, E_out, seed);	×
159	}
160
161	} // namespace openmc

openmc-dev / openmc / 17475377863

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