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

Use Discrete distribution for mesh sampling

Browse source 
This version builds the independent variable by
converting int to double and then converts back
when sampling.
This commit is contained in:
Paul P.H. Wilson 2023-03-16 17:34:44 -05:00
parent 2d20cd31f2
commit ee6975767d
2 changed files with 3 additions and 30 deletions
Download patch file Download diff file Expand all files Collapse all files

6
include/openmc/distribution_spatial.h
View file

@ -114,11 +114,7 @@ public:
private:
int32_t mesh_idx_ {C_NONE};
double total_strength_ {0.0};
// TODO: move to an independent class in the future that's similar
// to a discrete distribution without outcomes
std::vector<double> mesh_CDF_;
std::vector<double> mesh_strengths_;
UPtrDist elem_idx_; //!< Distribution of mesh element indices
};
//==============================================================================

27
src/distribution_spatial.cpp
View file

@ -212,10 +212,6 @@ MeshSpatial::MeshSpatial(pugi::xml_node node)
int32_t n_bins = this->n_sources();
std::vector<double> strengths(n_bins, 0.0);
mesh_CDF_.resize(n_bins + 1);
mesh_CDF_[0] = {0.0};
total_strength_ = 0.0;
// Create cdfs for sampling for an element over a mesh
// Volume scheme is weighted by the volume of each tet
// File scheme is weighted by an array given in the xml file
@ -228,28 +224,9 @@ MeshSpatial::MeshSpatial(pugi::xml_node node)
"not match the number of entities in mesh {} ({}).",
strengths.size(), mesh_id, n_bins));
}
mesh_strengths_ = std::move(strengths);
elem_idx_ = UPtrDist {new Discrete {strengths, n_bins}};
}
if (get_node_value_bool(node, "volume_normalized")) {
for (int i = 0; i < n_bins; i++) {
mesh_strengths_[i] *= mesh()->volume(i);
}
}
total_strength_ =
std::accumulate(mesh_strengths_.begin(), mesh_strengths_.end(), 0.0);
for (int i = 0; i < n_bins; i++) {
mesh_CDF_[i + 1] = mesh_CDF_[i] + mesh_strengths_[i] / total_strength_;
}
if (fabs(mesh_CDF_.back() - 1.0) > FP_COINCIDENT) {
fatal_error(
fmt::format("Mesh sampling CDF is incorrectly formed. Final value is: {}",
mesh_CDF_.back()));
}
mesh_CDF_.back() = 1.0;
}
Position MeshSpatial::sample(uint64_t* seed) const
@ -257,7 +234,7 @@ Position MeshSpatial::sample(uint64_t* seed) const
// Create random variable for sampling element from mesh
double eta = prn(seed);
// Sample over the CDF defined in initialization above
int32_t elem_idx = lower_bound_index(mesh_CDF_.begin(), mesh_CDF_.end(), eta);
int32_t elem_idx = elem_idx_->sample(eta);
return mesh()->sample(seed, elem_idx);
}