16008404600

Committed 01 Jul 2025 07:28PM UTC coverage: 85.223% (-0.1%) from 85.33%

Build # 16008404600

Build Type

Pull #3452

github

Committed by

web-flow

Commit Message

Merge 3621470de into eb74d497d

Pull Request Pull Request #3452: relaxed random ray constraints

Run Details

238 of 301 new or added lines in 8 files covered. (79.07%)

42 existing lines in 13 files now uncovered.

52846 of 62009 relevant lines covered (85.22%)

36646364.99 hits per line

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

74.06

/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 it) const
{
  return {dynamic_cast<FixedDistribution*>(x_.get())->sample(it),
    dynamic_cast<FixedDistribution*>(y_.get())->sample(it),
    dynamic_cast<FixedDistribution*>(z_.get())->sample(it)};
}

//==============================================================================
// 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);
  double z = z_->sample(seed);
  Position xi {r * cos(phi), r * sin(phi), z};
  return xi + origin_;
}
Position CylindricalIndependent::sample(vector<double>::iterator it) const
{
  double r = dynamic_cast<FixedDistribution*>(r_.get())->sample(it);
  double phi = dynamic_cast<FixedDistribution*>(phi_.get())->sample(it);
  double z = dynamic_cast<FixedDistribution*>(z_.get())->sample(it);
  Position xi {r * cos(phi), r * sin(phi), z};
  return xi + origin_;
}

//==============================================================================
// 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
  Position xi {std::sqrt(1 - cos_theta * cos_theta) * cos(phi),
    std::sqrt(1 - cos_theta * cos_theta) * sin(phi), cos_theta};
  return r * xi + origin_;
}
Position SphericalIndependent::sample(vector<double>::iterator it) const
{
  double r = dynamic_cast<FixedDistribution*>(r_.get())->sample(it);
  double cos_theta =
    dynamic_cast<FixedDistribution*>(cos_theta_.get())->sample(it);
  double phi = dynamic_cast<FixedDistribution*>(phi_.get())->sample(it);
  // sin(theta) by sin**2 + cos**2 = 1
  Position xi {std::sqrt(1 - cos_theta * cos_theta) * cos(phi),
    std::sqrt(1 - cos_theta * cos_theta) * sin(phi), cos_theta};
  return r * xi + origin_;
}

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

openmc-dev / openmc / 16008404600

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