17445237309

Committed 03 Sep 2025 08:25PM UTC coverage: 84.864% (-0.3%) from 85.209%

Build # 17445237309

Build Type

Pull #3460

github

Committed by

web-flow

Commit Message

Merge 5e7a2eae9 into 591856472

Pull Request Pull Request #3460: Source biasing capabilities

Run Details

343 of 638 new or added lines in 11 files covered. (53.76%)

4 existing lines in 3 files now uncovered.

53134 of 62611 relevant lines covered (84.86%)

38459716.64 hits per line

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

68.26

/src/distribution_spatial.cpp

#include "openmc/distribution_spatial.h"

#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/xml_interface.h"

namespace openmc {

//==============================================================================
// SpatialDistribution implementation
//==============================================================================

unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)
{
  // Check for type of spatial distribution and read
  std::string type;
  if (check_for_node(node, "type"))
    type = get_node_value(node, "type", true, true);
  if (type == "cartesian") {
    return UPtrSpace {new CartesianIndependent(node)};
  } else if (type == "cylindrical") {
    return UPtrSpace {new CylindricalIndependent(node)};
  } else if (type == "spherical") {
    return UPtrSpace {new SphericalIndependent(node)};
  } else if (type == "mesh") {
    return UPtrSpace {new MeshSpatial(node)};
  } else if (type == "cloud") {
    return UPtrSpace {new PointCloud(node)};
  } else if (type == "box") {
    return UPtrSpace {new SpatialBox(node)};
  } else if (type == "fission") {
    return UPtrSpace {new SpatialBox(node, true)};
  } else if (type == "point") {
    return UPtrSpace {new SpatialPoint(node)};
  } else {
    fatal_error(fmt::format(
      "Invalid spatial distribution for external source: {}", type));
  }
}

//==============================================================================
// CartesianIndependent implementation
//==============================================================================

CartesianIndependent::CartesianIndependent(pugi::xml_node node)
{
  // Read distribution for x coordinate
  if (check_for_node(node, "x")) {
    pugi::xml_node node_dist = node.child("x");
    x_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at x=0
    double x[] {0.0};
    double p[] {1.0};
    x_ = UPtrDist {new Discrete {x, p, 1}};
  }

  // Read distribution for y coordinate
  if (check_for_node(node, "y")) {
    pugi::xml_node node_dist = node.child("y");
    y_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at y=0
    double x[] {0.0};
    double p[] {1.0};
    y_ = UPtrDist {new Discrete {x, p, 1}};
  }

  // Read distribution for z coordinate
  if (check_for_node(node, "z")) {
    pugi::xml_node node_dist = node.child("z");
    z_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at z=0
    double x[] {0.0};
    double p[] {1.0};
    z_ = UPtrDist {new Discrete {x, p, 1}};
  }
}

std::pair<Position, double> CartesianIndependent::sample(uint64_t* seed) const
{
  auto [x_val, x_wgt] = x_->sample(seed);
  auto [y_val, y_wgt] = y_->sample(seed);
  auto [z_val, z_wgt] = z_->sample(seed);
  Position xi {x_val, y_val, z_val};
  return {xi, x_wgt * y_wgt * z_wgt};
}

//==============================================================================
// CylindricalIndependent implementation
//==============================================================================

CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
{
  // Read distribution for r-coordinate
  if (check_for_node(node, "r")) {
    pugi::xml_node node_dist = node.child("r");
    r_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at r=0
    double x[] {0.0};
    double p[] {1.0};
    r_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for phi-coordinate
  if (check_for_node(node, "phi")) {
    pugi::xml_node node_dist = node.child("phi");
    phi_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at phi=0
    double x[] {0.0};
    double p[] {1.0};
    phi_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for z-coordinate
  if (check_for_node(node, "z")) {
    pugi::xml_node node_dist = node.child("z");
    z_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at z=0
    double x[] {0.0};
    double p[] {1.0};
    z_ = make_unique<Discrete>(x, p, 1);
  }

  // Read cylinder center coordinates
  if (check_for_node(node, "origin")) {
    auto origin = get_node_array<double>(node, "origin");
    if (origin.size() == 3) {
      origin_ = origin;
    } else {
      fatal_error(
        "Origin for cylindrical source distribution must be length 3");
    }
  } else {
    // If no coordinates were specified, default to (0, 0, 0)
    origin_ = {0.0, 0.0, 0.0};
  }
}

std::pair<Position, double> CylindricalIndependent::sample(uint64_t* seed) const
{
  auto [r, r_wgt] = r_->sample(seed);
  auto [phi, phi_wgt] = phi_->sample(seed);
  auto [z, z_wgt] = z_->sample(seed);
  double x = r * cos(phi) + origin_.x;
  double y = r * sin(phi) + origin_.y;
  z += origin_.z;
  Position xi {x, y, z};
  return {xi, r_wgt * phi_wgt * z_wgt};
}

//==============================================================================
// SphericalIndependent implementation
//==============================================================================

SphericalIndependent::SphericalIndependent(pugi::xml_node node)
{
  // Read distribution for r-coordinate
  if (check_for_node(node, "r")) {
    pugi::xml_node node_dist = node.child("r");
    r_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at r=0
    double x[] {0.0};
    double p[] {1.0};
    r_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for cos_theta-coordinate
  if (check_for_node(node, "cos_theta")) {
    pugi::xml_node node_dist = node.child("cos_theta");
    cos_theta_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at
    // cos_theta=0
    double x[] {0.0};
    double p[] {1.0};
    cos_theta_ = make_unique<Discrete>(x, p, 1);
  }

  // Read distribution for phi-coordinate
  if (check_for_node(node, "phi")) {
    pugi::xml_node node_dist = node.child("phi");
    phi_ = distribution_from_xml(node_dist);
  } else {
    // If no distribution was specified, default to a single point at phi=0
    double x[] {0.0};
    double p[] {1.0};
    phi_ = make_unique<Discrete>(x, p, 1);
  }

  // Read sphere center coordinates
  if (check_for_node(node, "origin")) {
    auto origin = get_node_array<double>(node, "origin");
    if (origin.size() == 3) {
      origin_ = origin;
    } else {
      fatal_error("Origin for spherical source distribution must be length 3");
    }
  } else {
    // If no coordinates were specified, default to (0, 0, 0)
    origin_ = {0.0, 0.0, 0.0};
  }
}

std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const
{
  auto [r, r_wgt] = r_->sample(seed);
  auto [cos_theta, cos_theta_wgt] = cos_theta_->sample(seed);
  auto [phi, phi_wgt] = phi_->sample(seed);
  // sin(theta) by sin**2 + cos**2 = 1
  double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;
  double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
  double z = r * cos_theta + origin_.z;
  Position xi {x, y, z};
  return {xi, r_wgt * cos_theta_wgt * phi_wgt};
}

//==============================================================================
// MeshSpatial implementation
//==============================================================================

MeshSpatial::MeshSpatial(pugi::xml_node node)
{

  if (get_node_value(node, "type", true, true) != "mesh") {
    fatal_error(fmt::format(
      "Incorrect spatial type '{}' for a MeshSpatial distribution"));
  }

  // No in-tet distributions implemented, could include distributions for the
  // barycentric coords Read in unstructured mesh from mesh_id value
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id"));
  // Get pointer to spatial distribution
  mesh_idx_ = model::mesh_map.at(mesh_id);

  const auto mesh_ptr = model::meshes.at(mesh_idx_).get();

  check_element_types();

  size_t n_bins = this->n_sources();
  std::vector<double> strengths(n_bins, 1.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
  if (check_for_node(node, "strengths")) {
    strengths = get_node_array<double>(node, "strengths");
    if (strengths.size() != n_bins) {
      fatal_error(
        fmt::format("Number of entries in the source strengths array {} does "
                    "not match the number of entities in mesh {} ({}).",
          strengths.size(), mesh_id, n_bins));
    }
  }

  if (get_node_value_bool(node, "volume_normalized")) {
    for (int i = 0; i < n_bins; i++) {
      strengths[i] *= this->mesh()->volume(i);
    }
  }

  elem_idx_dist_.assign(strengths);

  if (check_for_node(node, "bias")) {
    pugi::xml_node bias_node = node.child("bias");

    if (check_for_node(bias_node, "strengths")) {
      std::vector<double> bias_strengths(n_bins, 1.0);
      bias_strengths = get_node_array<double>(node, "strengths");

      if (bias_strengths.size() != n_bins) {
        fatal_error(
          fmt::format("Number of entries in the bias strengths array {} does "
                      "not match the number of entities in mesh {} ({}).",
            bias_strengths.size(), mesh_id, n_bins));
      }

      if (get_node_value_bool(node, "volume_normalized")) {
        for (int i = 0; i < n_bins; i++) {
          bias_strengths[i] *= this->mesh()->volume(i);
        }
      }

      span<const double> b {bias_strengths};
      elem_idx_dist_.apply_bias(b);
    } else {
      fatal_error(fmt::format(
        "Bias node for mesh {} found without strengths array.", mesh_id));
    }
  }
}

MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)
  : mesh_idx_(mesh_idx)
{
  check_element_types();
  elem_idx_dist_.assign(strengths);
}

void MeshSpatial::check_element_types() const
{
  const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());
  if (umesh_ptr) {
    // ensure that the unstructured mesh contains only linear tets
    for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) {
      if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) {
        fatal_error(
          "Mesh specified for source must contain only linear tetrahedra.");
      }
    }
  }
}

int32_t MeshSpatial::sample_element_index(uint64_t* seed) const
{
  return elem_idx_dist_.sample(seed);
}

std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const
{
  // Sample the CDF defined in initialization above
  int32_t elem_idx = this->sample_element_index(seed);
  return {elem_idx, mesh()->sample_element(elem_idx, seed)};
}

std::pair<Position, double> MeshSpatial::sample(uint64_t* seed) const
{
  auto [elem_idx, u] = this->sample_mesh(seed);
  return {u, elem_idx_dist_.weight()[elem_idx]};
}

//==============================================================================
// PointCloud implementation
//==============================================================================

PointCloud::PointCloud(pugi::xml_node node)
{
  if (check_for_node(node, "coords")) {
    point_cloud_ = get_node_position_array(node, "coords");
  } else {
    fatal_error("No coordinates were provided for the PointCloud "
                "spatial distribution");
  }

  std::vector<double> strengths;

  if (check_for_node(node, "strengths"))
    strengths = get_node_array<double>(node, "strengths");
  else
    strengths.resize(point_cloud_.size(), 1.0);

  if (strengths.size() != point_cloud_.size()) {
    fatal_error(
      fmt::format("Number of entries for the strengths array {} does "
                  "not match the number of spatial points provided {}.",
        strengths.size(), point_cloud_.size()));
  }

  point_idx_dist_.assign(strengths);

  if (check_for_node(node, "bias")) {
    pugi::xml_node bias_node = node.child("bias");

    if (check_for_node(bias_node, "strengths")) {
      std::vector<double> bias_strengths(point_cloud_.size(), 1.0);
      bias_strengths = get_node_array<double>(node, "strengths");

      if (bias_strengths.size() != point_cloud_.size()) {
        fatal_error(
          fmt::format("Number of entries in the bias strengths array {} does "
                      "not match the number of spatial points provided {}.",
            bias_strengths.size(), point_cloud_.size()));
      }

      span<const double> b {bias_strengths};
      point_idx_dist_.apply_bias(b);
    } else {
      fatal_error(
        fmt::format("Bias node for PointCloud found without strengths array."));
    }
  }
}

PointCloud::PointCloud(
  std::vector<Position> point_cloud, span<const double> strengths)
{
  point_cloud_.assign(point_cloud.begin(), point_cloud.end());
  point_idx_dist_.assign(strengths);
}

std::pair<Position, double> PointCloud::sample(uint64_t* seed) const
{
  int32_t index = point_idx_dist_.sample(seed);
  return {point_cloud_[index], point_idx_dist_.weight()[index]};
}

//==============================================================================
// SpatialBox implementation
//==============================================================================

SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
  : only_fissionable_ {fission}
{
  // Read lower-right/upper-left coordinates
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 6)
    openmc::fatal_error("Box/fission spatial source must have six "
                        "parameters specified.");

  lower_left_ = Position {params[0], params[1], params[2]};
  upper_right_ = Position {params[3], params[4], params[5]};
}

std::pair<Position, double> SpatialBox::sample(uint64_t* seed) const
{
  Position xi {prn(seed), prn(seed), prn(seed)};
  // Biasing not implemented--use CartesianIndependent instead
  return {lower_left_ + xi * (upper_right_ - lower_left_), 1.0};
}

//==============================================================================
// SpatialPoint implementation
//==============================================================================

SpatialPoint::SpatialPoint(pugi::xml_node node)
{
  // Read location of point source
  auto params = get_node_array<double>(node, "parameters");
  if (params.size() != 3)
    openmc::fatal_error("Point spatial source must have three "
                        "parameters specified.");

  // Set position
  r_ = Position {params.data()};
}

std::pair<Position, double> SpatialPoint::sample(uint64_t* seed) const
{
  return {r_, 1.0};
}

} // namespace openmc

1	#include "openmc/distribution_spatial.h"
2
3	#include "openmc/error.h"
4	#include "openmc/mesh.h"
5	#include "openmc/random_lcg.h"
6	#include "openmc/search.h"
7	#include "openmc/xml_interface.h"
8
9	namespace openmc {
10
11	//==============================================================================
12	// SpatialDistribution implementation
13	//==============================================================================
14
15	unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)	6,802✔
16	{
17	// Check for type of spatial distribution and read
18	std::string type;	6,802✔
19	if (check_for_node(node, "type"))	6,802✔
20	type = get_node_value(node, "type", true, true);	6,802✔
21	if (type == "cartesian") {	6,802✔
22	return UPtrSpace {new CartesianIndependent(node)};	16✔
23	} else if (type == "cylindrical") {	6,786✔
24	return UPtrSpace {new CylindricalIndependent(node)};	48✔
25	} else if (type == "spherical") {	6,738✔
26	return UPtrSpace {new SphericalIndependent(node)};	63✔
27	} else if (type == "mesh") {	6,675✔
28	return UPtrSpace {new MeshSpatial(node)};	51✔
29	} else if (type == "cloud") {	6,624✔
30	return UPtrSpace {new PointCloud(node)};	14✔
31	} else if (type == "box") {	6,610✔
32	return UPtrSpace {new SpatialBox(node)};	3,735✔
33	} else if (type == "fission") {	2,875✔
34	return UPtrSpace {new SpatialBox(node, true)};	26✔
35	} else if (type == "point") {	2,849✔
36	return UPtrSpace {new SpatialPoint(node)};	2,849✔
37	} else {
38	fatal_error(fmt::format(	×
39	"Invalid spatial distribution for external source: {}", type));
40	}
41	}	6,801✔
42
43	//==============================================================================
44	// CartesianIndependent implementation
45	//==============================================================================
46
47	CartesianIndependent::CartesianIndependent(pugi::xml_node node)	16✔
48	{
49	// Read distribution for x coordinate
50	if (check_for_node(node, "x")) {	16✔
51	pugi::xml_node node_dist = node.child("x");	16✔
52	x_ = distribution_from_xml(node_dist);	16✔
53	} else {
54	// If no distribution was specified, default to a single point at x=0
55	double x[] {0.0};	×
56	double p[] {1.0};	×
57	x_ = UPtrDist {new Discrete {x, p, 1}};	×
58	}
59
60	// Read distribution for y coordinate
61	if (check_for_node(node, "y")) {	16✔
62	pugi::xml_node node_dist = node.child("y");	16✔
63	y_ = distribution_from_xml(node_dist);	16✔
64	} else {
65	// If no distribution was specified, default to a single point at y=0
66	double x[] {0.0};	×
67	double p[] {1.0};	×
68	y_ = UPtrDist {new Discrete {x, p, 1}};	×
69	}
70
71	// Read distribution for z coordinate
72	if (check_for_node(node, "z")) {	16✔
73	pugi::xml_node node_dist = node.child("z");	16✔
74	z_ = distribution_from_xml(node_dist);	16✔
75	} else {
76	// If no distribution was specified, default to a single point at z=0
77	double x[] {0.0};	×
78	double p[] {1.0};	×
79	z_ = UPtrDist {new Discrete {x, p, 1}};	×
80	}
81	}	16✔
82
83	std::pair<Position, double> CartesianIndependent::sample(uint64_t* seed) const	3,190✔
84	{
85	auto [x_val, x_wgt] = x_->sample(seed);	3,190✔
86	auto [y_val, y_wgt] = y_->sample(seed);	3,190✔
87	auto [z_val, z_wgt] = z_->sample(seed);	3,190✔
88	Position xi {x_val, y_val, z_val};	3,190✔
89	return {xi, x_wgt * y_wgt * z_wgt};	3,190✔
90	}
91
92	//==============================================================================
93	// CylindricalIndependent implementation
94	//==============================================================================
95
96	CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)	48✔
97	{
98	// Read distribution for r-coordinate
99	if (check_for_node(node, "r")) {	48✔
100	pugi::xml_node node_dist = node.child("r");	48✔
101	r_ = distribution_from_xml(node_dist);	48✔
102	} else {
103	// If no distribution was specified, default to a single point at r=0
104	double x[] {0.0};	×
105	double p[] {1.0};	×
106	r_ = make_unique<Discrete>(x, p, 1);	×
107	}
108
109	// Read distribution for phi-coordinate
110	if (check_for_node(node, "phi")) {	48✔
111	pugi::xml_node node_dist = node.child("phi");	48✔
112	phi_ = distribution_from_xml(node_dist);	48✔
113	} else {
114	// If no distribution was specified, default to a single point at phi=0
115	double x[] {0.0};	×
116	double p[] {1.0};	×
117	phi_ = make_unique<Discrete>(x, p, 1);	×
118	}
119
120	// Read distribution for z-coordinate
121	if (check_for_node(node, "z")) {	48✔
122	pugi::xml_node node_dist = node.child("z");	48✔
123	z_ = distribution_from_xml(node_dist);	48✔
124	} else {
125	// If no distribution was specified, default to a single point at z=0
126	double x[] {0.0};	×
127	double p[] {1.0};	×
128	z_ = make_unique<Discrete>(x, p, 1);	×
129	}
130
131	// Read cylinder center coordinates
132	if (check_for_node(node, "origin")) {	48✔
133	auto origin = get_node_array<double>(node, "origin");	48✔
134	if (origin.size() == 3) {	48✔
135	origin_ = origin;	48✔
136	} else {
137	fatal_error(	×
138	"Origin for cylindrical source distribution must be length 3");
139	}
140	} else {	48✔
141	// If no coordinates were specified, default to (0, 0, 0)
142	origin_ = {0.0, 0.0, 0.0};	×
143	}
144	}	48✔
145
146	std::pair<Position, double> CylindricalIndependent::sample(uint64_t* seed) const	3,542✔
147	{
148	auto [r, r_wgt] = r_->sample(seed);	3,542✔
149	auto [phi, phi_wgt] = phi_->sample(seed);	3,542✔
150	auto [z, z_wgt] = z_->sample(seed);	3,542✔
151	double x = r * cos(phi) + origin_.x;	3,542✔
152	double y = r * sin(phi) + origin_.y;	3,542✔
153	z += origin_.z;	3,542✔
154	Position xi {x, y, z};	3,542✔
155	return {xi, r_wgt * phi_wgt * z_wgt};	3,542✔
156	}
157
158	//==============================================================================
159	// SphericalIndependent implementation
160	//==============================================================================
161
162	SphericalIndependent::SphericalIndependent(pugi::xml_node node)	63✔
163	{
164	// Read distribution for r-coordinate
165	if (check_for_node(node, "r")) {	63✔
166	pugi::xml_node node_dist = node.child("r");	63✔
167	r_ = distribution_from_xml(node_dist);	63✔
168	} else {
169	// If no distribution was specified, default to a single point at r=0
170	double x[] {0.0};	×
171	double p[] {1.0};	×
172	r_ = make_unique<Discrete>(x, p, 1);	×
173	}
174
175	// Read distribution for cos_theta-coordinate
176	if (check_for_node(node, "cos_theta")) {	63✔
177	pugi::xml_node node_dist = node.child("cos_theta");	63✔
178	cos_theta_ = distribution_from_xml(node_dist);	63✔
179	} else {
180	// If no distribution was specified, default to a single point at
181	// cos_theta=0
182	double x[] {0.0};	×
183	double p[] {1.0};	×
184	cos_theta_ = make_unique<Discrete>(x, p, 1);	×
185	}
186
187	// Read distribution for phi-coordinate
188	if (check_for_node(node, "phi")) {	63✔
189	pugi::xml_node node_dist = node.child("phi");	63✔
190	phi_ = distribution_from_xml(node_dist);	63✔
191	} else {
192	// If no distribution was specified, default to a single point at phi=0
193	double x[] {0.0};	×
194	double p[] {1.0};	×
195	phi_ = make_unique<Discrete>(x, p, 1);	×
196	}
197
198	// Read sphere center coordinates
199	if (check_for_node(node, "origin")) {	63✔
200	auto origin = get_node_array<double>(node, "origin");	63✔
201	if (origin.size() == 3) {	63✔
202	origin_ = origin;	63✔
203	} else {
204	fatal_error("Origin for spherical source distribution must be length 3");	×
205	}
206	} else {	63✔
207	// If no coordinates were specified, default to (0, 0, 0)
208	origin_ = {0.0, 0.0, 0.0};	×
209	}
210	}	63✔
211
212	std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const	182,288✔
213	{
214	auto [r, r_wgt] = r_->sample(seed);	182,288✔
215	auto [cos_theta, cos_theta_wgt] = cos_theta_->sample(seed);	182,288✔
216	auto [phi, phi_wgt] = phi_->sample(seed);	182,288✔
217	// sin(theta) by sin2 + cos2 = 1
218	double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;	182,288✔
219	double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;	182,288✔
220	double z = r * cos_theta + origin_.z;	182,288✔
221	Position xi {x, y, z};	182,288✔
222	return {xi, r_wgt * cos_theta_wgt * phi_wgt};	182,288✔
223	}
224
225	//==============================================================================
226	// MeshSpatial implementation
227	//==============================================================================
228
229	MeshSpatial::MeshSpatial(pugi::xml_node node)	51✔
230	{
231
232	if (get_node_value(node, "type", true, true) != "mesh") {	51✔
233	fatal_error(fmt::format(	×
234	"Incorrect spatial type '{}' for a MeshSpatial distribution"));
235	}
236
237	// No in-tet distributions implemented, could include distributions for the
238	// barycentric coords Read in unstructured mesh from mesh_id value
239	int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id"));	51✔
240	// Get pointer to spatial distribution
241	mesh_idx_ = model::mesh_map.at(mesh_id);	51✔
242
243	const auto mesh_ptr = model::meshes.at(mesh_idx_).get();	51✔
244
245	check_element_types();	51✔
246
247	size_t n_bins = this->n_sources();	51✔
248	std::vector<double> strengths(n_bins, 1.0);	51✔
249
250	// Create cdfs for sampling for an element over a mesh
251	// Volume scheme is weighted by the volume of each tet
252	// File scheme is weighted by an array given in the xml file
253	if (check_for_node(node, "strengths")) {	51✔
254	strengths = get_node_array<double>(node, "strengths");	47✔
255	if (strengths.size() != n_bins) {	47✔
256	fatal_error(	1✔
257	fmt::format("Number of entries in the source strengths array {} does "	1✔
258	"not match the number of entities in mesh {} ({}).",
259	strengths.size(), mesh_id, n_bins));	1✔
260	}
261	}
262
263	if (get_node_value_bool(node, "volume_normalized")) {	50✔
264	for (int i = 0; i < n_bins; i++) {	49,726✔
265	strengths[i] *= this->mesh()->volume(i);	49,680✔
266	}
267	}
268
269	elem_idx_dist_.assign(strengths);	50✔
270
271	if (check_for_node(node, "bias")) {	50✔
NEW 272	pugi::xml_node bias_node = node.child("bias");	×
273
NEW 274	if (check_for_node(bias_node, "strengths")) {	×
NEW 275	std::vector<double> bias_strengths(n_bins, 1.0);	×
NEW 276	bias_strengths = get_node_array<double>(node, "strengths");	×
277
NEW 278	if (bias_strengths.size() != n_bins) {	×
NEW 279	fatal_error(	×
NEW 280	fmt::format("Number of entries in the bias strengths array {} does "	×
281	"not match the number of entities in mesh {} ({}).",
NEW 282	bias_strengths.size(), mesh_id, n_bins));	×
283	}
284
NEW 285	if (get_node_value_bool(node, "volume_normalized")) {	×
NEW 286	for (int i = 0; i < n_bins; i++) {	×
NEW 287	bias_strengths[i] *= this->mesh()->volume(i);	×
288	}
289	}
290
NEW 291	span<const double> b {bias_strengths};	×
NEW 292	elem_idx_dist_.apply_bias(b);	×
NEW 293	} else {	×
NEW 294	fatal_error(fmt::format(	×
295	"Bias node for mesh {} found without strengths array.", mesh_id));
296	}
297	}
298	}	50✔
299
300	MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)	208✔
301	: mesh_idx_(mesh_idx)	208✔
302	{
303	check_element_types();	208✔
304	elem_idx_dist_.assign(strengths);	208✔
305	}	208✔
306
307	void MeshSpatial::check_element_types() const	259✔
308	{
309	const auto umesh_ptr = dynamic_cast<const UnstructuredMesh*>(this->mesh());	259✔
310	if (umesh_ptr) {	259✔
311	// ensure that the unstructured mesh contains only linear tets
312	for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) {	156,013✔
313	if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) {	156,000✔
314	fatal_error(	×
315	"Mesh specified for source must contain only linear tetrahedra.");
316	}
317	}
318	}
319	}	259✔
320
321	int32_t MeshSpatial::sample_element_index(uint64_t* seed) const	2,679,368✔
322	{
323	return elem_idx_dist_.sample(seed);	2,679,368✔
324	}
325
326	std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const	802,020✔
327	{
328	// Sample the CDF defined in initialization above
329	int32_t elem_idx = this->sample_element_index(seed);	802,020✔
330	return {elem_idx, mesh()->sample_element(elem_idx, seed)};	802,020✔
331	}
332
333	std::pair<Position, double> MeshSpatial::sample(uint64_t* seed) const	802,020✔
334	{
335	auto [elem_idx, u] = this->sample_mesh(seed);	802,020✔
336	return {u, elem_idx_dist_.weight()[elem_idx]};	1,604,040✔
337	}
338
339	//==============================================================================
340	// PointCloud implementation
341	//==============================================================================
342
343	PointCloud::PointCloud(pugi::xml_node node)	14✔
344	{
345	if (check_for_node(node, "coords")) {	14✔
346	point_cloud_ = get_node_position_array(node, "coords");	14✔
347	} else {
348	fatal_error("No coordinates were provided for the PointCloud "	×
349	"spatial distribution");
350	}
351
352	std::vector<double> strengths;	14✔
353
354	if (check_for_node(node, "strengths"))	14✔
355	strengths = get_node_array<double>(node, "strengths");	14✔
356	else
357	strengths.resize(point_cloud_.size(), 1.0);	×
358
359	if (strengths.size() != point_cloud_.size()) {	14✔
360	fatal_error(	×
361	fmt::format("Number of entries for the strengths array {} does "	×
362	"not match the number of spatial points provided {}.",
363	strengths.size(), point_cloud_.size()));	×
364	}
365
366	point_idx_dist_.assign(strengths);	14✔
367
368	if (check_for_node(node, "bias")) {	14✔
NEW 369	pugi::xml_node bias_node = node.child("bias");	×
370
NEW 371	if (check_for_node(bias_node, "strengths")) {	×
NEW 372	std::vector<double> bias_strengths(point_cloud_.size(), 1.0);	×
NEW 373	bias_strengths = get_node_array<double>(node, "strengths");	×
374
NEW 375	if (bias_strengths.size() != point_cloud_.size()) {	×
NEW 376	fatal_error(	×
NEW 377	fmt::format("Number of entries in the bias strengths array {} does "	×
378	"not match the number of spatial points provided {}.",
NEW 379	bias_strengths.size(), point_cloud_.size()));	×
380	}
381
NEW 382	span<const double> b {bias_strengths};	×
NEW 383	point_idx_dist_.apply_bias(b);	×
NEW 384	} else {	×
NEW 385	fatal_error(	×
NEW 386	fmt::format("Bias node for PointCloud found without strengths array."));	×
387	}
388	}
389	}	14✔
390
391	PointCloud::PointCloud(	×
392	std::vector<Position> point_cloud, span<const double> strengths)	×
393	{
394	point_cloud_.assign(point_cloud.begin(), point_cloud.end());	×
395	point_idx_dist_.assign(strengths);	×
396	}
397
398	std::pair<Position, double> PointCloud::sample(uint64_t* seed) const	700,000✔
399	{
400	int32_t index = point_idx_dist_.sample(seed);	700,000✔
401	return {point_cloud_[index], point_idx_dist_.weight()[index]};	700,000✔
402	}
403
404	//==============================================================================
405	// SpatialBox implementation
406	//==============================================================================
407
408	SpatialBox::SpatialBox(pugi::xml_node node, bool fission)	3,761✔
409	: only_fissionable_ {fission}	3,761✔
410	{
411	// Read lower-right/upper-left coordinates
412	auto params = get_node_array<double>(node, "parameters");	3,761✔
413	if (params.size() != 6)	3,761✔
414	openmc::fatal_error("Box/fission spatial source must have six "	×
415	"parameters specified.");
416
417	lower_left_ = Position {params[0], params[1], params[2]};	3,761✔
418	upper_right_ = Position {params[3], params[4], params[5]};	3,761✔
419	}	3,761✔
420
421	std::pair<Position, double> SpatialBox::sample(uint64_t* seed) const	5,236,921✔
422	{
423	Position xi {prn(seed), prn(seed), prn(seed)};	5,236,921✔
424	// Biasing not implemented--use CartesianIndependent instead
425	return {lower_left_ + xi * (upper_right_ - lower_left_), 1.0};	5,236,921✔
426	}
427
428	//==============================================================================
429	// SpatialPoint implementation
430	//==============================================================================
431
432	SpatialPoint::SpatialPoint(pugi::xml_node node)	2,849✔
433	{
434	// Read location of point source
435	auto params = get_node_array<double>(node, "parameters");	2,849✔
436	if (params.size() != 3)	2,849✔
437	openmc::fatal_error("Point spatial source must have three "	×
438	"parameters specified.");
439
440	// Set position
441	r_ = Position {params.data()};	2,849✔
442	}	2,849✔
443
444	std::pair<Position, double> SpatialPoint::sample(uint64_t* seed) const	26,462,596✔
445	{
446	return {r_, 1.0};	26,462,596✔
447	}
448
449	} // namespace openmc

openmc-dev / openmc / 17445237309

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