15900894277

Committed 26 Jun 2025 11:42AM UTC coverage: 85.214% (-0.04%) from 85.252%

Build # 15900894277

Build Type

Pull #3452

github

Committed by

web-flow

Commit Message

Merge a81fd95f9 into 5c1021446

Pull Request Pull Request #3452: relaxed random ray constraints

Run Details

237 of 306 new or added lines in 8 files covered. (77.45%)

2 existing lines in 1 file now uncovered.

52791 of 61951 relevant lines covered (85.21%)

36510182.52 hits per line

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

72.13

/src/distribution_spatial.cpp

#include "openmc/distribution_spatial.h"

#include <cmath>

#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/mesh.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/vector.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}};
  }
  if ((x_->dims() > 0) && (y_->dims() > 0) && (z_->dims() > 0))
    dims_ = x_->dims() + y_->dims() + z_->dims();
  volume_ = 1.0;
  auto sup = x_->support();
  volume_ *= sup.second - sup.first;
  sup = y_->support();
  volume_ *= sup.second - sup.first;
  sup = z_->support();
  volume_ *= sup.second - sup.first;
}

Position CartesianIndependent::sample(uint64_t* seed) const
{
  return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};
}
Position CartesianIndependent::sample(vector<double>::iterator x) const
{
  const auto* _x = dynamic_cast<FixedDistribution*>(x_.get());
  const auto* _y = dynamic_cast<FixedDistribution*>(y_.get());
  const auto* _z = dynamic_cast<FixedDistribution*>(z_.get());
  return {_x->sample(x), _y->sample(x), _z->sample(x)};
}

//==============================================================================
// 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};
  }
  if ((r_->dims() > 0) && (phi_->dims() > 0) && (z_->dims() > 0))
    dims_ = r_->dims() + phi_->dims() + z_->dims();
  volume_ = 1.0;
  auto sup = r_->support();
  volume_ *= 0.5 * (std::pow(sup.second, 2) - std::pow(sup.first, 2));
  sup = phi_->support();
  volume_ *= sup.second - sup.first;
  sup = z_->support();
  volume_ *= sup.second - sup.first;
}

Position CylindricalIndependent::sample(uint64_t* seed) const
{
  double r = r_->sample(seed);
  double phi = phi_->sample(seed);
  return {r * cos(phi) + origin_.x, r * sin(phi) + origin_.y,
    z_->sample(seed) + origin_.z};
}
Position CylindricalIndependent::sample(vector<double>::iterator x) const
{
  const auto* _r = dynamic_cast<FixedDistribution*>(r_.get());
  const auto* _phi = dynamic_cast<FixedDistribution*>(phi_.get());
  const auto* _z = dynamic_cast<FixedDistribution*>(z_.get());
  double r = _r->sample(x);
  double phi = _phi->sample(x);
  return {r * cos(phi) + origin_.x, r * sin(phi) + origin_.y,
    _z->sample(x) + origin_.z};
}

//==============================================================================
// 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};
  }
  if ((r_->dims() > 0) && (cos_theta_->dims() > 0) && (phi_->dims() > 0))
    dims_ = r_->dims() + cos_theta_->dims() + phi_->dims();
  volume_ = 1.0;
  auto sup = r_->support();
  volume_ *= 1.0 / 3.0 * (std::pow(sup.second, 3) - std::pow(sup.first, 3));
  sup = cos_theta_->support();
  volume_ *= sup.second - sup.first;
  sup = phi_->support();
  volume_ *= sup.second - sup.first;
}

Position SphericalIndependent::sample(uint64_t* seed) const
{
  double r = r_->sample(seed);
  double cos_theta = cos_theta_->sample(seed);
  double phi = phi_->sample(seed);
  // sin(theta) by sin**2 + cos**2 = 1
  return {r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x,
    r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y,
    r * cos_theta + origin_.z};
}
Position SphericalIndependent::sample(vector<double>::iterator x) const
{
  const auto* _r = dynamic_cast<FixedDistribution*>(r_.get());
  const auto* _cos_theta = dynamic_cast<FixedDistribution*>(cos_theta_.get());
  const auto* _phi = dynamic_cast<FixedDistribution*>(phi_.get());
  double r = _r->sample(x);
  double cos_theta = _cos_theta->sample(x);
  double phi = _phi->sample(x);
  // sin(theta) by sin**2 + cos**2 = 1
  return {r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x,
    r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y,
    r * cos_theta + origin_.z};
}

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

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)};
}

Position MeshSpatial::sample(uint64_t* seed) const
{
  return this->sample_mesh(seed).second;
}

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

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);
}

Position PointCloud::sample(uint64_t* seed) const
{
  int32_t index = point_idx_dist_.sample(seed);
  return point_cloud_[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]};
  dims_ = 3;
  if (!only_fissionable_)
    volume_ =
      ((upper_right_[0] - lower_left_[0]) * (upper_right_[1] - lower_left_[1]) *
        (upper_right_[2] - lower_left_[2]));
}

Position SpatialBox::sample(vector<double>::iterator x) const
{
  Position xi {*(x++), *(x++), *(x++)};
  return lower_left_ + xi * (upper_right_ - lower_left_);
}

//==============================================================================
// 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()};
}

Position SpatialPoint::sample(uint64_t* seed) const
{
  return r_;
}

} // namespace openmc

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

openmc-dev / openmc / 15900894277

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