8745330834

Committed 18 Apr 2024 10:10PM UTC coverage: 84.629% (-0.01%) from 84.643%

Build # 8745330834

Build Type

push

github

Committed by

web-flow

Commit Message

Random Ray Transport (#2823)

Co-authored-by: Gavin Ridley <gavin.keith.ridley@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Run Details

577 of 700 new or added lines in 16 files covered. (82.43%)

3 existing lines in 2 files now uncovered.

48402 of 57193 relevant lines covered (84.63%)

34065627.92 hits per line

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

85.06

/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/mcpl_interface.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/string_utils.h"
#include "openmc/xml_interface.h"

namespace openmc {

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

namespace model {

vector<unique_ptr<Source>> external_sources;
}

//==============================================================================
// Source create implementation
//==============================================================================

unique_ptr<Source> Source::create(pugi::xml_node node)
{
  // if the source type is present, use it to determine the type
  // of object to create
  if (check_for_node(node, "type")) {
    std::string source_type = get_node_value(node, "type");
    if (source_type == "independent") {
      return make_unique<IndependentSource>(node);
    } else if (source_type == "file") {
      return make_unique<FileSource>(node);
    } else if (source_type == "compiled") {
      return make_unique<CompiledSourceWrapper>(node);
    } else if (source_type == "mesh") {
      return make_unique<MeshSource>(node);
    } else {
      fatal_error(fmt::format("Invalid source type '{}' found.", source_type));
    }
  } else {
    // support legacy source format
    if (check_for_node(node, "file")) {
      return make_unique<FileSource>(node);
    } else if (check_for_node(node, "library")) {
      return make_unique<CompiledSourceWrapper>(node);
    } else {
      return make_unique<IndependentSource>(node);
    }
  }
}

//==============================================================================
// 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")) {
      space_ = SpatialDistribution::create(node.child("space"));
    } 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")) {
      angle_ = UnitSphereDistribution::create(node.child("angle"));
    } 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);

  // Sample energy and time for neutron and photon sources
  if (settings::solver_type != SolverType::RANDOM_RAY) {
    // 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(pugi::xml_node node)
{
  auto path = get_node_value(node, "file", false, true);
  if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) {
    sites_ = mcpl_source_sites(path);
  } else {
    this->load_sites_from_file(path);
  }
}

FileSource::FileSource(const std::string& path)
{
  load_sites_from_file(path);
}

void FileSource::load_sites_from_file(const 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(pugi::xml_node node)
{
  // Get shared library path and parameters
  auto path = get_node_value(node, "library", false, true);
  std::string parameters;
  if (check_for_node(node, "parameters")) {
    parameters = get_node_value(node, "parameters", false, true);
  }
  setup(path, parameters);
}

void CompiledSourceWrapper::setup(
  const std::string& path, const 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
}

//==============================================================================
// MeshSource implementation
//==============================================================================

MeshSource::MeshSource(pugi::xml_node node)
{
  int32_t mesh_id = stoi(get_node_value(node, "mesh"));
  int32_t mesh_idx = model::mesh_map.at(mesh_id);
  const auto& mesh = model::meshes[mesh_idx];

  std::vector<double> strengths;
  // read all source distributions and populate strengths vector for MeshSpatial
  // object
  for (auto source_node : node.children("source")) {
    sources_.emplace_back(Source::create(source_node));
    strengths.push_back(sources_.back()->strength());
  }

  // the number of source distributions should either be one or equal to the
  // number of mesh elements
  if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) {
    fatal_error(fmt::format("Incorrect number of source distributions ({}) for "
                            "mesh source with {} elements.",
      sources_.size(), mesh->n_bins()));
  }

  space_ = std::make_unique<MeshSpatial>(mesh_idx, strengths);
}

SourceSite MeshSource::sample(uint64_t* seed) const
{
  // sample location and element from mesh
  auto mesh_location = space_->sample_mesh(seed);

  // Sample source for the chosen element
  int32_t element = mesh_location.first;
  SourceSite site = source(element)->sample(seed);

  // Replace spatial position with the one already sampled
  site.r = mesh_location.second;

  return site;
}

//==============================================================================
// 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/mcpl_interface.h"
26	#include "openmc/memory.h"
27	#include "openmc/message_passing.h"
28	#include "openmc/mgxs_interface.h"
29	#include "openmc/nuclide.h"
30	#include "openmc/random_lcg.h"
31	#include "openmc/search.h"
32	#include "openmc/settings.h"
33	#include "openmc/simulation.h"
34	#include "openmc/state_point.h"
35	#include "openmc/string_utils.h"
36	#include "openmc/xml_interface.h"
37
38	namespace openmc {
39
40	//==============================================================================
41	// Global variables
42	//==============================================================================
43
44	namespace model {
45
46	vector<unique_ptr<Source>> external_sources;
47	}
48
49	//==============================================================================
50	// Source create implementation
51	//==============================================================================
52
53	unique_ptr<Source> Source::create(pugi::xml_node node)	44,272✔
54	{
55	// if the source type is present, use it to determine the type
56	// of object to create
57	if (check_for_node(node, "type")) {	44,272✔
58	std::string source_type = get_node_value(node, "type");	43,004✔
59	if (source_type == "independent") {	43,004✔
60	return make_unique<IndependentSource>(node);	42,685✔
61	} else if (source_type == "file") {	319✔
62	return make_unique<FileSource>(node);	40✔
63	} else if (source_type == "compiled") {	279✔
64	return make_unique<CompiledSourceWrapper>(node);	38✔
65	} else if (source_type == "mesh") {	241✔
66	return make_unique<MeshSource>(node);	241✔
67	} else {
68	fatal_error(fmt::format("Invalid source type '{}' found.", source_type));	×
69	}
70	} else {	42,991✔
71	// support legacy source format
72	if (check_for_node(node, "file")) {	1,268✔
73	return make_unique<FileSource>(node);	38✔
74	} else if (check_for_node(node, "library")) {	1,230✔
75	return make_unique<CompiledSourceWrapper>(node);	×
76	} else {
77	return make_unique<IndependentSource>(node);	1,230✔
78	}
79	}
80	}
81
82	//==============================================================================
83	// IndependentSource implementation
84	//==============================================================================
85
86	IndependentSource::IndependentSource(	1,790✔
87	UPtrSpace space, UPtrAngle angle, UPtrDist energy, UPtrDist time)	1,790✔
88	: space_ {std::move(space)}, angle_ {std::move(angle)},	1,790✔
89	energy_ {std::move(energy)}, time_ {std::move(time)}	3,580✔
90	{}	1,790✔
91
92	IndependentSource::IndependentSource(pugi::xml_node node)	43,915✔
93	{
94	// Check for particle type
95	if (check_for_node(node, "particle")) {	43,915✔
96	auto temp_str = get_node_value(node, "particle", true, true);	42,685✔
97	if (temp_str == "neutron") {	42,685✔
98	particle_ = ParticleType::neutron;	42,572✔
99	} else if (temp_str == "photon") {	113✔
100	particle_ = ParticleType::photon;	113✔
101	settings::photon_transport = true;	113✔
102	} else {
103	fatal_error(std::string("Unknown source particle type: ") + temp_str);	×
104	}
105	}	42,685✔
106
107	// Check for source strength
108	if (check_for_node(node, "strength")) {	43,915✔
109	strength_ = std::stod(get_node_value(node, "strength"));	43,117✔
110	}
111
112	// Check for external source file
113	if (check_for_node(node, "file")) {	43,915✔
114
115	} else {
116
117	// Spatial distribution for external source
118	if (check_for_node(node, "space")) {	43,915✔
119	space_ = SpatialDistribution::create(node.child("space"));	6,029✔
120	} else {
121	// If no spatial distribution specified, make it a point source
122	space_ = UPtrSpace {new SpatialPoint()};	37,886✔
123	}
124
125	// Determine external source angular distribution
126	if (check_for_node(node, "angle")) {	43,914✔
127	angle_ = UnitSphereDistribution::create(node.child("angle"));	3,949✔
128	} else {
129	angle_ = UPtrAngle {new Isotropic()};	39,965✔
130	}
131
132	// Determine external source energy distribution
133	if (check_for_node(node, "energy")) {	43,914✔
134	pugi::xml_node node_dist = node.child("energy");	4,473✔
135	energy_ = distribution_from_xml(node_dist);	4,473✔
136	} else {
137	// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
138	energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};	39,441✔
139	}
140
141	// Determine external source time distribution
142	if (check_for_node(node, "time")) {	43,914✔
143	pugi::xml_node node_dist = node.child("time");	52✔
144	time_ = distribution_from_xml(node_dist);	52✔
145	} else {
146	// Default to a Constant time T=0
147	double T[] {0.0};	43,862✔
148	double p[] {1.0};	43,862✔
149	time_ = UPtrDist {new Discrete {T, p, 1}};	43,862✔
150	}
151
152	// Check for domains to reject from
153	if (check_for_node(node, "domain_type")) {	43,914✔
154	std::string domain_type = get_node_value(node, "domain_type");	14✔
155	if (domain_type == "cell") {	14✔
156	domain_type_ = DomainType::CELL;	14✔
157	} else if (domain_type == "material") {	×
158	domain_type_ = DomainType::MATERIAL;	×
159	} else if (domain_type == "universe") {	×
160	domain_type_ = DomainType::UNIVERSE;	×
161	} else {
162	fatal_error(std::string(	×
163	"Unrecognized domain type for source rejection: " + domain_type));
164	}
165
166	auto ids = get_node_array<int>(node, "domain_ids");	14✔
167	domain_ids_.insert(ids.begin(), ids.end());	14✔
168	}	14✔
169	}
170	}	43,914✔
171
172	SourceSite IndependentSource::sample(uint64_t* seed) const	17,379,775✔
173	{
174	SourceSite site;	17,379,775✔
175	site.particle = particle_;	17,379,775✔
176
177	// Repeat sampling source location until a good site has been found
178	bool found = false;	17,379,775✔
179	int n_reject = 0;	17,379,775✔
180	static int n_accept = 0;
181
182	while (!found) {	37,259,560✔
183	// Set particle type
184	Particle p;	19,879,785✔
185	p.type() = particle_;	19,879,785✔
186	p.u() = {0.0, 0.0, 1.0};	19,879,785✔
187
188	// Sample spatial distribution
189	p.r() = space_->sample(seed);	19,879,785✔
190
191	// Now search to see if location exists in geometry
192	found = exhaustive_find_cell(p);	19,879,785✔
193
194	// Check if spatial site is in fissionable material
195	if (found) {	19,879,785✔
196	auto space_box = dynamic_cast<SpatialBox*>(space_.get());	19,433,661✔
197	if (space_box) {	19,433,661✔
198	if (space_box->only_fissionable()) {	4,853,880✔
199	// Determine material
200	auto mat_index = p.material();	1,192,482✔
201	if (mat_index == MATERIAL_VOID) {	1,192,482✔
202	found = false;	×
203	} else {
204	found = model::materials[mat_index]->fissionable();	1,192,482✔
205	}
206	}
207	}
208
209	// Rejection based on cells/materials/universes
210	if (!domain_ids_.empty()) {	19,433,661✔
211	found = false;	1,741,404✔
212	if (domain_type_ == DomainType::MATERIAL) {	1,741,404✔
213	auto mat_index = p.material();	×
214	if (mat_index != MATERIAL_VOID) {	×
215	found = contains(domain_ids_, model::materials[mat_index]->id());	×
216	}
217	} else {
218	for (int i = 0; i < p.n_coord(); i++) {	3,468,808✔
219	auto id = (domain_type_ == DomainType::CELL)	1,741,404✔
220	? model::cells[p.coord(i).cell]->id_	1,741,404✔
221	: model::universes[p.coord(i).universe]->id_;	×
222	if ((found = contains(domain_ids_, id)))	1,741,404✔
223	break;	14,000✔
224	}
225	}
226	}
227	}
228
229	// Check for rejection
230	if (!found) {	19,879,785✔
231	++n_reject;	2,500,010✔
232	if (n_reject >= EXTSRC_REJECT_THRESHOLD &&	2,500,010✔
233	static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
234	fatal_error("More than 95% of external source sites sampled were "	×
235	"rejected. Please check your external source's spatial "
236	"definition.");
237	}
238	}
239
240	site.r = p.r();	19,879,785✔
241	}	19,879,785✔
242
243	// Sample angle
244	site.u = angle_->sample(seed);	17,379,775✔
245
246	// Sample energy and time for neutron and photon sources
247	if (settings::solver_type != SolverType::RANDOM_RAY) {	17,379,775✔
248	// Check for monoenergetic source above maximum particle energy
249	auto p = static_cast<int>(particle_);	17,351,775✔
250	auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());	17,351,775✔
251	if (energy_ptr) {	17,351,775✔
252	auto energies = xt::adapt(energy_ptr->x());	6,378,180✔
253	if (xt::any(energies > data::energy_max[p])) {	6,378,180✔
NEW 254	fatal_error("Source energy above range of energies of at least "	×
255	"one cross section table");
256	}
257	}	6,378,180✔
258
259	while (true) {
260	// Sample energy spectrum
261	site.E = energy_->sample(seed);	17,351,775✔
262
263	// Resample if energy falls above maximum particle energy
264	if (site.E < data::energy_max[p])	17,351,775✔
265	break;	17,351,775✔
266
NEW 267	n_reject++;	×
NEW 268	if (n_reject >= EXTSRC_REJECT_THRESHOLD &&	×
269	static_cast<double>(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) {
NEW 270	fatal_error(	×
271	"More than 95% of external source sites sampled were "
272	"rejected. Please check your external source energy spectrum "
273	"definition.");
274	}
275	}
276
277	// Sample particle creation time
278	site.time = time_->sample(seed);	17,351,775✔
279	}
280
281	// Increment number of accepted samples
282	++n_accept;	17,379,775✔
283
284	return site;	17,379,775✔
285	}
286
287	//==============================================================================
288	// FileSource implementation
289	//==============================================================================
290	FileSource::FileSource(pugi::xml_node node)	78✔
291	{
292	auto path = get_node_value(node, "file", false, true);	78✔
293	if (ends_with(path, ".mcpl") \|\| ends_with(path, ".mcpl.gz")) {	78✔
294	sites_ = mcpl_source_sites(path);	19✔
295	} else {
296	this->load_sites_from_file(path);	59✔
297	}
298	}	66✔
299
300	FileSource::FileSource(const std::string& path)	19✔
301	{
302	load_sites_from_file(path);	19✔
303	}	19✔
304
305	void FileSource::load_sites_from_file(const std::string& path)	78✔
306	{
307	// Check if source file exists
308	if (!file_exists(path)) {	78✔
309	fatal_error(fmt::format("Source file '{}' does not exist.", path));	×
310	}
311
312	// Read the source from a binary file instead of sampling from some
313	// assumed source distribution
314	write_message(6, "Reading source file from {}...", path);	78✔
315
316	// Open the binary file
317	hid_t file_id = file_open(path, 'r', true);	78✔
318
319	// Check to make sure this is a source file
320	std::string filetype;	78✔
321	read_attribute(file_id, "filetype", filetype);	78✔
322	if (filetype != "source" && filetype != "statepoint") {	78✔
323	fatal_error("Specified starting source file not a source file type.");	×
324	}
325
326	// Read in the source particles
327	read_source_bank(file_id, sites_, false);	78✔
328
329	// Close file
330	file_close(file_id);	66✔
331	}	66✔
332
333	SourceSite FileSource::sample(uint64_t* seed) const	168,560✔
334	{
335	size_t i_site = sites_.size() * prn(seed);	168,560✔
336	return sites_[i_site];	168,560✔
337	}
338
339	//==============================================================================
340	// CompiledSourceWrapper implementation
341	//==============================================================================
342	CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node)	38✔
343	{
344	// Get shared library path and parameters
345	auto path = get_node_value(node, "library", false, true);	38✔
346	std::string parameters;	38✔
347	if (check_for_node(node, "parameters")) {	38✔
348	parameters = get_node_value(node, "parameters", false, true);	19✔
349	}
350	setup(path, parameters);	38✔
351	}	38✔
352
353	void CompiledSourceWrapper::setup(	38✔
354	const std::string& path, const std::string& parameters)
355	{
356	#ifdef HAS_DYNAMIC_LINKING
357	// Open the library
358	shared_library_ = dlopen(path.c_str(), RTLD_LAZY);	38✔
359	if (!shared_library_) {	38✔
360	fatal_error("Couldn't open source library " + path);	×
361	}
362
363	// reset errors
364	dlerror();	38✔
365
366	// get the function to create the custom source from the library
367	auto create_compiled_source = reinterpret_cast<create_compiled_source_t*>(
368	dlsym(shared_library_, "openmc_create_source"));	38✔
369
370	// check for any dlsym errors
371	auto dlsym_error = dlerror();	38✔
372	if (dlsym_error) {	38✔
373	std::string error_msg = fmt::format(
374	"Couldn't open the openmc_create_source symbol: {}", dlsym_error);	×
375	dlclose(shared_library_);	×
376	fatal_error(error_msg);	×
377	}	×
378
379	// create a pointer to an instance of the custom source
380	compiled_source_ = create_compiled_source(parameters);	38✔
381
382	#else
383	fatal_error("Custom source libraries have not yet been implemented for "
384	"non-POSIX systems");
385	#endif
386	}	38✔
387
388	CompiledSourceWrapper::~CompiledSourceWrapper()	76✔
389	{
390	// Make sure custom source is cleared before closing shared library
391	if (compiled_source_.get())	38✔
392	compiled_source_.reset();	38✔
393
394	#ifdef HAS_DYNAMIC_LINKING
395	dlclose(shared_library_);	38✔
396	#else
397	fatal_error("Custom source libraries have not yet been implemented for "
398	"non-POSIX systems");
399	#endif
400	}	76✔
401		38✔
402	//==============================================================================
403	// MeshSource implementation
404	//==============================================================================
405
406	MeshSource::MeshSource(pugi::xml_node node)
407	{
408	int32_t mesh_id = stoi(get_node_value(node, "mesh"));
409	int32_t mesh_idx = model::mesh_map.at(mesh_id);
410	const auto& mesh = model::meshes[mesh_idx];
411
412	std::vector<double> strengths;
413	// read all source distributions and populate strengths vector for MeshSpatial	38✔
414	// object	38✔
415	for (auto source_node : node.children("source")) {
416	sources_.emplace_back(Source::create(source_node));
417	strengths.push_back(sources_.back()->strength());	38✔
418	}	38✔
419
420	// the number of source distributions should either be one or equal to the
421	// number of mesh elements	38✔
422	if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) {
423	fatal_error(fmt::format("Incorrect number of source distributions ({}) for "
424	"mesh source with {} elements.",
425	sources_.size(), mesh->n_bins()));
426	}	38✔
427
428	space_ = std::make_unique<MeshSpatial>(mesh_idx, strengths);
429	}
430
431	SourceSite MeshSource::sample(uint64_t* seed) const
432	{	241✔
433	// sample location and element from mesh
434	auto mesh_location = space_->sample_mesh(seed);	241✔
435		241✔
436	// Sample source for the chosen element	241✔
437	int32_t element = mesh_location.first;
438	SourceSite site = source(element)->sample(seed);	241✔
439
440	// Replace spatial position with the one already sampled
441	site.r = mesh_location.second;	38,047✔
442		37,806✔
443	return site;	37,806✔
444	}
445
446	//==============================================================================
447	// Non-member functions
448	//==============================================================================	241✔
449		×
450	void initialize_source()
451	{	×
452	write_message("Initializing source particles...", 5);
453
454	// Generation source sites from specified distribution in user input	241✔
455	#pragma omp parallel for	241✔
456	for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
457	// initialize random number seed	224,170✔
458	int64_t id = simulation::total_gen * settings::n_particles +
459	simulation::work_index[mpi::rank] + i + 1;
460	uint64_t seed = init_seed(id, STREAM_SOURCE);	224,170✔
461
462	// sample external source distribution
463	simulation::source_bank[i] = sample_external_source(&seed);	224,170✔
464	}	224,170✔
465
466	// Write out initial source
467	if (settings::write_initial_source) {	224,170✔
468	write_message("Writing out initial source...", 5);
469	std::string filename = settings::path_output + "initial_source.h5";	448,340✔
470	hid_t file_id = file_open(filename, 'w', true);
471	write_source_bank(file_id, simulation::source_bank, simulation::work_index);
472	file_close(file_id);
473	}
474	}
475
476	SourceSite sample_external_source(uint64_t* seed)	3,901✔
477	{
478	// Determine total source strength	3,901✔
479	double total_strength = 0.0;
480	for (auto& s : model::external_sources)
481	total_strength += s->strength();
482		1,648,743✔
483	// Sample from among multiple source distributions
484	int i = 0;	3,293,400✔
485	if (model::external_sources.size() > 1) {	1,646,700✔
486	double xi = prn(seed) * total_strength;	1,646,700✔
487	double c = 0.0;
488	for (; i < model::external_sources.size(); ++i) {
489	c += model::external_sources[i]->strength();	1,646,700✔
490	if (xi < c)
491	break;
492	}
493	}	3,901✔
494		×
495	// Sample source site from i-th source distribution	×
496	SourceSite site {model::external_sources[i]->sample(seed)};	×
497		×
498	// If running in MG, convert site.E to group	×
499	if (!settings::run_CE) {
500	site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),	3,901✔
501	data::mg.rev_energy_bins_.end(), site.E);
502	site.E = data::mg.num_energy_groups_ - site.E - 1.;	17,800,335✔
503	}
504
505	return site;	17,800,335✔
506	}	41,130,670✔
507		23,330,335✔
508	void free_memory_source()
509	{
510	model::external_sources.clear();	17,800,335✔
511	}	17,800,335✔
512		182,000✔
513	//==============================================================================	182,000✔
514	// C API	2,922,528✔
515	//==============================================================================	2,922,528✔
516		2,922,528✔
517	extern "C" int openmc_sample_external_source(	182,000✔
518	size_t n, uint64_t* seed, void* sites)
519	{
520	if (!sites \|\| !seed) {
521	set_errmsg("Received null pointer.");
522	return OPENMC_E_INVALID_ARGUMENT;	17,800,335✔
523	}
524
525	if (model::external_sources.empty()) {	17,800,335✔
526	set_errmsg("No external sources have been defined.");	817,600✔
527	return OPENMC_E_OUT_OF_BOUNDS;
528	}	817,600✔
529
530	auto sites_array = static_cast<SourceSite*>(sites);
531	for (size_t i = 0; i < n; ++i) {	17,800,335✔
532	sites_array[i] = sample_external_source(seed);
533	}
534	return 0;	6,602✔
535	}
536		6,602✔
537	} // namespace openmc	6,602✔

openmc-dev / openmc / 8745330834

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