6726911467

Committed 02 Nov 2023 12:56AM UTC coverage: 84.679% (-0.003%) from 84.682%

Build # 6726911467

Build Type

push

github

Committed by

web-flow

Commit Message

Fix loop over particle coordinate levels for source domain rejection (#2751)

Run Details

2 of 3 new or added lines in 1 file covered. (66.67%)

1 existing line in 1 file now uncovered.

47130 of 55657 relevant lines covered (84.68%)

85932242.34 hits per line

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

84.4

/src/source.cpp

#include "openmc/source.h"

#if defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#define HAS_DYNAMIC_LINKING
#endif

#include <algorithm> // for move

#ifdef HAS_DYNAMIC_LINKING
#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
#endif

#include "xtensor/xadapt.hpp"
#include <fmt/core.h>

#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/container_util.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/memory.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/state_point.h"
#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================
// Global variables
//==============================================================================

namespace model {

vector<unique_ptr<Source>> external_sources;
}

//==============================================================================
// IndependentSource implementation
//==============================================================================

IndependentSource::IndependentSource(
  UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time)
  : space_ {std::move(space)}, angle_ {std::move(angle)},
    energy_ {std::move(energy)}, time_ {std::move(time)}
{}

IndependentSource::IndependentSource(pugi::xml_node node)
{
  // Check for particle type
  if (check_for_node(node, "particle")) {
    auto temp_str = get_node_value(node, "particle", true, true);
    if (temp_str == "neutron") {
      particle_ = ParticleType::neutron;
    } else if (temp_str == "photon") {
      particle_ = ParticleType::photon;
      settings::photon_transport = true;
    } else {
      fatal_error(std::string("Unknown source particle type: ") + temp_str);
    }
  }

  // Check for source strength
  if (check_for_node(node, "strength")) {
    strength_ = std::stod(get_node_value(node, "strength"));
  }

  // Check for external source file
  if (check_for_node(node, "file")) {

  } else {

    // Spatial distribution for external source
    if (check_for_node(node, "space")) {
      // Get pointer to spatial distribution
      pugi::xml_node node_space = node.child("space");

      // Check for type of spatial distribution and read
      std::string type;
      if (check_for_node(node_space, "type"))
        type = get_node_value(node_space, "type", true, true);
      if (type == "cartesian") {
        space_ = UPtrSpace {new CartesianIndependent(node_space)};
      } else if (type == "cylindrical") {
        space_ = UPtrSpace {new CylindricalIndependent(node_space)};
      } else if (type == "spherical") {
        space_ = UPtrSpace {new SphericalIndependent(node_space)};
      } else if (type == "mesh") {
        space_ = UPtrSpace {new MeshSpatial(node_space)};
      } else if (type == "box") {
        space_ = UPtrSpace {new SpatialBox(node_space)};
      } else if (type == "fission") {
        space_ = UPtrSpace {new SpatialBox(node_space, true)};
      } else if (type == "point") {
        space_ = UPtrSpace {new SpatialPoint(node_space)};
      } else {
        fatal_error(fmt::format(
          "Invalid spatial distribution for external source: {}", type));
      }

    } else {
      // If no spatial distribution specified, make it a point source
      space_ = UPtrSpace {new SpatialPoint()};
    }

    // Determine external source angular distribution
    if (check_for_node(node, "angle")) {
      // Get pointer to angular distribution
      pugi::xml_node node_angle = node.child("angle");

      // Check for type of angular distribution
      std::string type;
      if (check_for_node(node_angle, "type"))
        type = get_node_value(node_angle, "type", true, true);
      if (type == "isotropic") {
        angle_ = UPtrAngle {new Isotropic()};
      } else if (type == "monodirectional") {
        angle_ = UPtrAngle {new Monodirectional(node_angle)};
      } else if (type == "mu-phi") {
        angle_ = UPtrAngle {new PolarAzimuthal(node_angle)};
      } else {
        fatal_error(fmt::format(
          "Invalid angular distribution for external source: {}", type));
      }

    } else {
      angle_ = UPtrAngle {new Isotropic()};
    }

    // Determine external source energy distribution
    if (check_for_node(node, "energy")) {
      pugi::xml_node node_dist = node.child("energy");
      energy_ = distribution_from_xml(node_dist);
    } else {
      // Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
      energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};
    }

    // Determine external source time distribution
    if (check_for_node(node, "time")) {
      pugi::xml_node node_dist = node.child("time");
      time_ = distribution_from_xml(node_dist);
    } else {
      // Default to a Constant time T=0
      double T[] {0.0};
      double p[] {1.0};
      time_ = UPtrDist {new Discrete {T, p, 1}};
    }

    // Check for domains to reject from
    if (check_for_node(node, "domain_type")) {
      std::string domain_type = get_node_value(node, "domain_type");
      if (domain_type == "cell") {
        domain_type_ = DomainType::CELL;
      } else if (domain_type == "material") {
        domain_type_ = DomainType::MATERIAL;
      } else if (domain_type == "universe") {
        domain_type_ = DomainType::UNIVERSE;
      } else {
        fatal_error(std::string(
          "Unrecognized domain type for source rejection: " + domain_type));
      }

      auto ids = get_node_array<int>(node, "domain_ids");
      domain_ids_.insert(ids.begin(), ids.end());
    }
  }
}

SourceSite IndependentSource::sample(uint64_t* seed) const
{
  SourceSite site;
  site.particle = particle_;

  // Repeat sampling source location until a good site has been found
  bool found = false;
  int n_reject = 0;
  static int n_accept = 0;

  while (!found) {
    // Set particle type
    Particle p;
    p.type() = particle_;
    p.u() = {0.0, 0.0, 1.0};

    // Sample spatial distribution
    p.r() = space_->sample(seed);

    // Now search to see if location exists in geometry
    found = exhaustive_find_cell(p);

    // Check if spatial site is in fissionable material
    if (found) {
      auto space_box = dynamic_cast<SpatialBox*>(space_.get());
      if (space_box) {
        if (space_box->only_fissionable()) {
          // Determine material
          auto mat_index = p.material();
          if (mat_index == MATERIAL_VOID) {
            found = false;
          } else {
            found = model::materials[mat_index]->fissionable_;
          }
        }
      }

      // Rejection based on cells/materials/universes
      if (!domain_ids_.empty()) {
        found = false;
        if (domain_type_ == DomainType::MATERIAL) {
          auto mat_index = p.material();
          if (mat_index != MATERIAL_VOID) {
            found = contains(domain_ids_, model::materials[mat_index]->id());
          }
        } else {
          for (int i = 0; i < p.n_coord(); i++) {
            auto id = (domain_type_ == DomainType::CELL)
                        ? model::cells[p.coord(i).cell]->id_
                        : model::universes[p.coord(i).universe]->id_;
            if ((found = contains(domain_ids_, id)))
              break;
          }
        }
      }
    }

    // Check for rejection
    if (!found) {
      ++n_reject;
      if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
          static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
        fatal_error("More than 95% of external source sites sampled were "
                    "rejected. Please check your external source's spatial "
                    "definition.");
      }
    }

    site.r = p.r();
  }

  // Sample angle
  site.u = angle_->sample(seed);

  // Check for monoenergetic source above maximum particle energy
  auto p = static_cast<int>(particle_);
  auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
  if (energy_ptr) {
    auto energies = xt::adapt(energy_ptr->x());
    if (xt::any(energies > data::energy_max[p])) {
      fatal_error("Source energy above range of energies of at least "
                  "one cross section table");
    }
  }

  while (true) {
    // Sample energy spectrum
    site.E = energy_->sample(seed);

    // Resample if energy falls above maximum particle energy
    if (site.E < data::energy_max[p])
      break;

    n_reject++;
    if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
        static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
      fatal_error("More than 95% of external source sites sampled were "
                  "rejected. Please check your external source energy spectrum "
                  "definition.");
    }
  }

  // Sample particle creation time
  site.time = time_->sample(seed);

  // Increment number of accepted samples
  ++n_accept;

  return site;
}

//==============================================================================
// FileSource implementation
//==============================================================================

FileSource::FileSource(std::string path)
{
  // Check if source file exists
  if (!file_exists(path)) {
    fatal_error(fmt::format("Source file '{}' does not exist.", path));
  }

  // Read the source from a binary file instead of sampling from some
  // assumed source distribution
  write_message(6, "Reading source file from {}...", path);

  // Open the binary file
  hid_t file_id = file_open(path, 'r', true);

  // Check to make sure this is a source file
  std::string filetype;
  read_attribute(file_id, "filetype", filetype);
  if (filetype != "source" && filetype != "statepoint") {
    fatal_error("Specified starting source file not a source file type.");
  }

  // Read in the source particles
  read_source_bank(file_id, sites_, false);

  // Close file
  file_close(file_id);
}

SourceSite FileSource::sample(uint64_t* seed) const
{
  size_t i_site = sites_.size() * prn(seed);
  return sites_[i_site];
}

//==============================================================================
// CompiledSourceWrapper implementation
//==============================================================================

CompiledSourceWrapper::CompiledSourceWrapper(
  std::string path, std::string parameters)
{
#ifdef HAS_DYNAMIC_LINKING
  // Open the library
  shared_library_ = dlopen(path.c_str(), RTLD_LAZY);
  if (!shared_library_) {
    fatal_error("Couldn't open source library " + path);
  }

  // reset errors
  dlerror();

  // get the function to create the custom source from the library
  auto create_compiled_source = reinterpret_cast<create_compiled_source_t*>(
    dlsym(shared_library_, "openmc_create_source"));

  // check for any dlsym errors
  auto dlsym_error = dlerror();
  if (dlsym_error) {
    std::string error_msg = fmt::format(
      "Couldn't open the openmc_create_source symbol: {}", dlsym_error);
    dlclose(shared_library_);
    fatal_error(error_msg);
  }

  // create a pointer to an instance of the custom source
  compiled_source_ = create_compiled_source(parameters);

#else
  fatal_error("Custom source libraries have not yet been implemented for "
              "non-POSIX systems");
#endif
}

CompiledSourceWrapper::~CompiledSourceWrapper()
{
  // Make sure custom source is cleared before closing shared library
  if (compiled_source_.get())
    compiled_source_.reset();

#ifdef HAS_DYNAMIC_LINKING
  dlclose(shared_library_);
#else
  fatal_error("Custom source libraries have not yet been implemented for "
              "non-POSIX systems");
#endif
}

//==============================================================================
// Non-member functions
//==============================================================================

void initialize_source()
{
  write_message("Initializing source particles...", 5);

// Generation source sites from specified distribution in user input
#pragma omp parallel for
  for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
    // initialize random number seed
    int64_t id = simulation::total_gen * settings::n_particles +
                 simulation::work_index[mpi::rank] + i + 1;
    uint64_t seed = init_seed(id, STREAM_SOURCE);

    // sample external source distribution
    simulation::source_bank[i] = sample_external_source(&seed);
  }

  // Write out initial source
  if (settings::write_initial_source) {
    write_message("Writing out initial source...", 5);
    std::string filename = settings::path_output + "initial_source.h5";
    hid_t file_id = file_open(filename, 'w', true);
    write_source_bank(file_id, simulation::source_bank, simulation::work_index);
    file_close(file_id);
  }
}

SourceSite sample_external_source(uint64_t* seed)
{
  // Determine total source strength
  double total_strength = 0.0;
  for (auto& s : model::external_sources)
    total_strength += s->strength();

  // Sample from among multiple source distributions
  int i = 0;
  if (model::external_sources.size() > 1) {
    double xi = prn(seed) * total_strength;
    double c = 0.0;
    for (; i < model::external_sources.size(); ++i) {
      c += model::external_sources[i]->strength();
      if (xi < c)
        break;
    }
  }

  // Sample source site from i-th source distribution
  SourceSite site {model::external_sources[i]->sample(seed)};

  // If running in MG, convert site.E to group
  if (!settings::run_CE) {
    site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),
      data::mg.rev_energy_bins_.end(), site.E);
    site.E = data::mg.num_energy_groups_ - site.E - 1.;
  }

  return site;
}

void free_memory_source()
{
  model::external_sources.clear();
}

//==============================================================================
// C API
//==============================================================================

extern "C" int openmc_sample_external_source(
  size_t n, uint64_t* seed, void* sites)
{
  if (!sites || !seed) {
    set_errmsg("Received null pointer.");
    return OPENMC_E_INVALID_ARGUMENT;
  }

  if (model::external_sources.empty()) {
    set_errmsg("No external sources have been defined.");
    return OPENMC_E_OUT_OF_BOUNDS;
  }

  auto sites_array = static_cast<SourceSite*>(sites);
  for (size_t i = 0; i < n; ++i) {
    sites_array[i] = sample_external_source(seed);
  }
  return 0;
}

} // namespace openmc

1	#include "openmc/source.h"
2
3	#if defined(__unix__) \|\| (defined(__APPLE__) && defined(__MACH__))
4	#define HAS_DYNAMIC_LINKING
5	#endif
6
7	#include <algorithm> // for move
8
9	#ifdef HAS_DYNAMIC_LINKING
10	#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
11	#endif
12
13	#include "xtensor/xadapt.hpp"
14	#include <fmt/core.h>
15
16	#include "openmc/bank.h"
17	#include "openmc/capi.h"
18	#include "openmc/cell.h"
19	#include "openmc/container_util.h"
20	#include "openmc/error.h"
21	#include "openmc/file_utils.h"
22	#include "openmc/geometry.h"
23	#include "openmc/hdf5_interface.h"
24	#include "openmc/material.h"
25	#include "openmc/memory.h"
26	#include "openmc/message_passing.h"
27	#include "openmc/mgxs_interface.h"
28	#include "openmc/nuclide.h"
29	#include "openmc/random_lcg.h"
30	#include "openmc/search.h"
31	#include "openmc/settings.h"
32	#include "openmc/simulation.h"
33	#include "openmc/state_point.h"
34	#include "openmc/xml_interface.h"
35
36	namespace openmc {
37
38	//==============================================================================
39	// Global variables
40	//==============================================================================
41
42	namespace model {
43
44	vector<unique_ptr<Source>> external_sources;
45	}
46
47	//==============================================================================
48	// IndependentSource implementation
49	//==============================================================================
50
51	IndependentSource::IndependentSource(	1,769✔
52	UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time)	1,769✔
53	: space_ {std::move(space)}, angle_ {std::move(angle)},	1,769✔
54	energy_ {std::move(energy)}, time_ {std::move(time)}	3,538✔
55	{}	1,769✔
56
57	IndependentSource::IndependentSource(pugi::xml_node node)	6,931✔
58	{
59	// Check for particle type
60	if (check_for_node(node, "particle")) {	6,931✔
61	auto temp_str = get_node_value(node, "particle", true, true);	4,786✔
62	if (temp_str == "neutron") {	4,786✔
63	particle_ = ParticleType::neutron;	4,673✔
64	} else if (temp_str == "photon") {	113✔
65	particle_ = ParticleType::photon;	113✔
66	settings::photon_transport = true;	113✔
67	} else {
68	fatal_error(std::string("Unknown source particle type: ") + temp_str);	×
69	}
70	}	4,786✔
71
72	// Check for source strength
73	if (check_for_node(node, "strength")) {	6,931✔
74	strength_ = std::stod(get_node_value(node, "strength"));	6,133✔
75	}
76
77	// Check for external source file
78	if (check_for_node(node, "file")) {	6,931✔
79
80	} else {
81
82	// Spatial distribution for external source
83	if (check_for_node(node, "space")) {	6,931✔
84	// Get pointer to spatial distribution
85	pugi::xml_node node_space = node.child("space");	6,837✔
86
87	// Check for type of spatial distribution and read
88	std::string type;	6,837✔
89	if (check_for_node(node_space, "type"))	6,837✔
90	type = get_node_value(node_space, "type", true, true);	6,837✔
91	if (type == "cartesian") {	6,837✔
92	space_ = UPtrSpace {new CartesianIndependent(node_space)};	19✔
93	} else if (type == "cylindrical") {	6,818✔
94	space_ = UPtrSpace {new CylindricalIndependent(node_space)};	57✔
95	} else if (type == "spherical") {	6,761✔
96	space_ = UPtrSpace {new SphericalIndependent(node_space)};	69✔
97	} else if (type == "mesh") {	6,692✔
98	space_ = UPtrSpace {new MeshSpatial(node_space)};	7✔
99	} else if (type == "box") {	6,685✔
100	space_ = UPtrSpace {new SpatialBox(node_space)};	2,720✔
101	} else if (type == "fission") {	3,965✔
102	space_ = UPtrSpace {new SpatialBox(node_space, true)};	1,011✔
103	} else if (type == "point") {	2,954✔
104	space_ = UPtrSpace {new SpatialPoint(node_space)};	2,954✔
105	} else {
106	fatal_error(fmt::format(	×
107	"Invalid spatial distribution for external source: {}", type));
108	}
109
110	} else {	6,836✔
111	// If no spatial distribution specified, make it a point source
112	space_ = UPtrSpace {new SpatialPoint()};	94✔
113	}
114
115	// Determine external source angular distribution
116	if (check_for_node(node, "angle")) {	6,930✔
117	// Get pointer to angular distribution
118	pugi::xml_node node_angle = node.child("angle");	2,157✔
119
120	// Check for type of angular distribution
121	std::string type;	2,157✔
122	if (check_for_node(node_angle, "type"))	2,157✔
123	type = get_node_value(node_angle, "type", true, true);	2,157✔
124	if (type == "isotropic") {	2,157✔
125	angle_ = UPtrAngle {new Isotropic()};	301✔
126	} else if (type == "monodirectional") {	1,856✔
127	angle_ = UPtrAngle {new Monodirectional(node_angle)};	1,823✔
128	} else if (type == "mu-phi") {	33✔
129	angle_ = UPtrAngle {new PolarAzimuthal(node_angle)};	33✔
130	} else {
131	fatal_error(fmt::format(	×
132	"Invalid angular distribution for external source: {}", type));
133	}
134
135	} else {	2,157✔
136	angle_ = UPtrAngle {new Isotropic()};	4,773✔
137	}
138
139	// Determine external source energy distribution
140	if (check_for_node(node, "energy")) {	6,930✔
141	pugi::xml_node node_dist = node.child("energy");	2,681✔
142	energy_ = distribution_from_xml(node_dist);	2,681✔
143	} else {
144	// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
145	energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};	4,249✔
146	}
147
148	// Determine external source time distribution
149	if (check_for_node(node, "time")) {	6,930✔
150	pugi::xml_node node_dist = node.child("time");	52✔
151	time_ = distribution_from_xml(node_dist);	52✔
152	} else {
153	// Default to a Constant time T=0
154	double T[] {0.0};	6,878✔
155	double p[] {1.0};	6,878✔
156	time_ = UPtrDist {new Discrete {T, p, 1}};	6,878✔
157	}
158
159	// Check for domains to reject from
160	if (check_for_node(node, "domain_type")) {	6,930✔
161	std::string domain_type = get_node_value(node, "domain_type");	14✔
162	if (domain_type == "cell") {	14✔
163	domain_type_ = DomainType::CELL;	14✔
164	} else if (domain_type == "material") {	×
165	domain_type_ = DomainType::MATERIAL;	×
166	} else if (domain_type == "universe") {	×
167	domain_type_ = DomainType::UNIVERSE;	×
168	} else {
169	fatal_error(std::string(	×
170	"Unrecognized domain type for source rejection: " + domain_type));
171	}
172
173	auto ids = get_node_array<int>(node, "domain_ids");	14✔
174	domain_ids_.insert(ids.begin(), ids.end());	14✔
175	}	14✔
176	}
177	}	6,930✔
178
179	SourceSite IndependentSource::sample(uint64_t* seed) const	17,585,245✔
180	{
181	SourceSite site;	17,585,245✔
182	site.particle = particle_;	17,585,245✔
183
184	// Repeat sampling source location until a good site has been found
185	bool found = false;	17,585,245✔
186	int n_reject = 0;	17,585,245✔
187	static int n_accept = 0;
188
189	while (!found) {	37,447,570✔
190	// Set particle type
191	Particle p;	19,862,325✔
192	p.type() = particle_;	19,862,325✔
193	p.u() = {0.0, 0.0, 1.0};	19,862,325✔
194
195	// Sample spatial distribution
196	p.r() = space_->sample(seed);	19,862,325✔
197
198	// Now search to see if location exists in geometry
199	found = exhaustive_find_cell(p);	19,862,325✔
200
201	// Check if spatial site is in fissionable material
202	if (found) {	19,862,325✔
203	auto space_box = dynamic_cast<SpatialBox*>(space_.get());	19,818,617✔
204	if (space_box) {	19,818,617✔
205	if (space_box->only_fissionable()) {	5,384,266✔
206	// Determine material
207	auto mat_index = p.material();	1,492,654✔
208	if (mat_index == MATERIAL_VOID) {	1,492,654✔
209	found = false;	×
210	} else {
211	found = model::materials[mat_index]->fissionable_;	1,492,654✔
212	}
213	}
214	}
215
216	// Rejection based on cells/materials/universes
217	if (!domain_ids_.empty()) {	19,818,617✔
218	found = false;	1,715,518✔
219	if (domain_type_ == DomainType::MATERIAL) {	1,715,518✔
220	auto mat_index = p.material();	×
221	if (mat_index != MATERIAL_VOID) {	×
222	found = contains(domain_ids_, model::materials[mat_index]->id());	×
223	}
224	} else {
225	for (int i = 0; i < p.n_coord(); i++) {	3,417,036✔
226	auto id = (domain_type_ == DomainType::CELL)	1,715,518✔
227	? model::cells[p.coord(i).cell]->id_	1,715,518✔
NEW 228	: model::universes[p.coord(i).universe]->id_;	×
229	if ((found = contains(domain_ids_, id)))	1,715,518✔
230	break;	14,000✔
231	}
232	}
233	}
234	}
235
236	// Check for rejection
237	if (!found) {	19,862,325✔
238	++n_reject;	2,277,080✔
239	if (n_reject >= EXTSRC_REJECT_THRESHOLD &&	2,277,080✔
240	static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
241	fatal_error("More than 95% of external source sites sampled were "	×
242	"rejected. Please check your external source's spatial "
243	"definition.");
244	}
245	}
246
247	site.r = p.r();	19,862,325✔
248	}	19,862,325✔
249
250	// Sample angle
251	site.u = angle_->sample(seed);	17,585,245✔
252
253	// Check for monoenergetic source above maximum particle energy
254	auto p = static_cast<int>(particle_);	17,585,245✔
255	auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());	17,585,245✔
256	if (energy_ptr) {	17,585,245✔
257	auto energies = xt::adapt(energy_ptr->x());	6,154,180✔
258	if (xt::any(energies > data::energy_max[p])) {	6,154,180✔
259	fatal_error("Source energy above range of energies of at least "	×
260	"one cross section table");
261	}
262	}	6,154,180✔
263
264	while (true) {
265	// Sample energy spectrum
266	site.E = energy_->sample(seed);	17,585,245✔
267
268	// Resample if energy falls above maximum particle energy
269	if (site.E < data::energy_max[p])	17,585,245✔
270	break;	17,585,245✔
271
272	n_reject++;	×
273	if (n_reject >= EXTSRC_REJECT_THRESHOLD &&	×
274	static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
275	fatal_error("More than 95% of external source sites sampled were "	×
276	"rejected. Please check your external source energy spectrum "
277	"definition.");
278	}
279	}
280
281	// Sample particle creation time
282	site.time = time_->sample(seed);	17,585,245✔
283
284	// Increment number of accepted samples
285	++n_accept;	17,585,245✔
286
287	return site;	17,585,245✔
288	}
289
290	//==============================================================================
291	// FileSource implementation
292	//==============================================================================
293
294	FileSource::FileSource(std::string path)	64✔
295	{
296	// Check if source file exists
297	if (!file_exists(path)) {	64✔
298	fatal_error(fmt::format("Source file '{}' does not exist.", path));	×
299	}
300
301	// Read the source from a binary file instead of sampling from some
302	// assumed source distribution
303	write_message(6, "Reading source file from {}...", path);	64✔
304
305	// Open the binary file
306	hid_t file_id = file_open(path, 'r', true);	64✔
307
308	// Check to make sure this is a source file
309	std::string filetype;	64✔
310	read_attribute(file_id, "filetype", filetype);	64✔
311	if (filetype != "source" && filetype != "statepoint") {	64✔
312	fatal_error("Specified starting source file not a source file type.");	×
313	}
314
315	// Read in the source particles
316	read_source_bank(file_id, sites_, false);	64✔
317
318	// Close file
319	file_close(file_id);	52✔
320	}	52✔
321
322	SourceSite FileSource::sample(uint64_t* seed) const	168,420✔
323	{
324	size_t i_site = sites_.size() * prn(seed);	168,420✔
325	return sites_[i_site];	168,420✔
326	}
327
328	//==============================================================================
329	// CompiledSourceWrapper implementation
330	//==============================================================================
331
332	CompiledSourceWrapper::CompiledSourceWrapper(	38✔
333	std::string path, std::string parameters)	38✔
334	{
335	#ifdef HAS_DYNAMIC_LINKING
336	// Open the library
337	shared_library_ = dlopen(path.c_str(), RTLD_LAZY);	38✔
338	if (!shared_library_) {	38✔
339	fatal_error("Couldn't open source library " + path);	×
340	}
341
342	// reset errors
343	dlerror();	38✔
344
345	// get the function to create the custom source from the library
346	auto create_compiled_source = reinterpret_cast<create_compiled_source_t*>(
347	dlsym(shared_library_, "openmc_create_source"));	38✔
348
349	// check for any dlsym errors
350	auto dlsym_error = dlerror();	38✔
351	if (dlsym_error) {	38✔
352	std::string error_msg = fmt::format(
353	"Couldn't open the openmc_create_source symbol: {}", dlsym_error);	×
354	dlclose(shared_library_);	×
355	fatal_error(error_msg);	×
356	}	×
357
358	// create a pointer to an instance of the custom source
359	compiled_source_ = create_compiled_source(parameters);	38✔
360
361	#else
362	fatal_error("Custom source libraries have not yet been implemented for "
363	"non-POSIX systems");
364	#endif
365	}	38✔
366
367	CompiledSourceWrapper::~CompiledSourceWrapper()	76✔
368	{
369	// Make sure custom source is cleared before closing shared library
370	if (compiled_source_.get())	38✔
371	compiled_source_.reset();	38✔
372
373	#ifdef HAS_DYNAMIC_LINKING
374	dlclose(shared_library_);	38✔
375	#else
376	fatal_error("Custom source libraries have not yet been implemented for "
377	"non-POSIX systems");
378	#endif
379	}	76✔
380		38✔
381	//==============================================================================
382	// Non-member functions
383	//==============================================================================
384
385	void initialize_source()
386	{
387	write_message("Initializing source particles...", 5);
388
389	// Generation source sites from specified distribution in user input
390	#pragma omp parallel for
391	for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
392	// initialize random number seed	38✔
393	int64_t id = simulation::total_gen * settings::n_particles +	38✔
394	simulation::work_index[mpi::rank] + i + 1;
395	uint64_t seed = init_seed(id, STREAM_SOURCE);
396		38✔
397	// sample external source distribution	38✔
398	simulation::source_bank[i] = sample_external_source(&seed);
399	}
400		38✔
401	// Write out initial source
402	if (settings::write_initial_source) {
403	write_message("Writing out initial source...", 5);
404	std::string filename = settings::path_output + "initial_source.h5";
405	hid_t file_id = file_open(filename, 'w', true);	38✔
406	write_source_bank(file_id, simulation::source_bank, simulation::work_index);
407	file_close(file_id);
408	}
409	}
410
411	SourceSite sample_external_source(uint64_t* seed)	4,769✔
412	{
413	// Determine total source strength	4,769✔
414	double total_strength = 0.0;
415	for (auto& s : model::external_sources)
416	total_strength += s->strength();
417		1,905,167✔
418	// Sample from among multiple source distributions
419	int i = 0;	3,805,800✔
420	if (model::external_sources.size() > 1) {	1,902,900✔
421	double xi = prn(seed) * total_strength;	1,902,900✔
422	double c = 0.0;
423	for (; i < model::external_sources.size(); ++i) {
424	c += model::external_sources[i]->strength();	1,902,900✔
425	if (xi < c)
426	break;
427	}
428	}	4,769✔
429		×
430	// Sample source site from i-th source distribution	×
431	SourceSite site {model::external_sources[i]->sample(seed)};	×
432		×
433	// If running in MG, convert site.E to group	×
434	if (!settings::run_CE) {
435	site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),	4,769✔
436	data::mg.rev_energy_bins_.end(), site.E);
437	site.E = data::mg.num_energy_groups_ - site.E - 1.;	18,033,665✔
438	}
439
440	return site;	18,033,665✔
441	}	41,597,330✔
442		23,563,665✔
443	void free_memory_source()
444	{
445	model::external_sources.clear();	18,033,665✔
446	}	18,033,665✔
447		182,000✔
448	//==============================================================================	182,000✔
449	// C API	2,922,528✔
450	//==============================================================================	2,922,528✔
451		2,922,528✔
452	extern "C" int openmc_sample_external_source(	182,000✔
453	size_t n, uint64_t* seed, void* sites)
454	{
455	if (!sites \|\| !seed) {
456	set_errmsg("Received null pointer.");
457	return OPENMC_E_INVALID_ARGUMENT;	18,033,665✔
458	}
459
460	if (model::external_sources.empty()) {	18,033,665✔
461	set_errmsg("No external sources have been defined.");	817,600✔
462	return OPENMC_E_OUT_OF_BOUNDS;
463	}	817,600✔
464
465	auto sites_array = static_cast<SourceSite*>(sites);
466	for (size_t i = 0; i < n; ++i) {	18,033,665✔
467	sites_array[i] = sample_external_source(seed);
468	}
469	return 0;	7,171✔
470	}
471		7,171✔
472	} // namespace openmc	7,171✔

openmc-dev / openmc / 6726911467

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