21892714514

Committed 10 Feb 2026 11:58PM UTC coverage: 74.806% (-0.02%) from 74.828%

Build # 21892714514

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #929 from wdhawkins/iteration_log

Fixing iteration status messages

Coverage Stats

4 of 4 new or added lines in 1 file covered. (100.0%)

211 existing lines in 16 files now uncovered.

19638 of 26252 relevant lines covered (74.81%)

69654996.54 hits per line

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

82.61

/modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.cc

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/math/functions/function.h"
#include "framework/object_factory.h"
#include "framework/logging/log.h"
#include "framework/runtime.h"
#include <numeric>
#include <limits>

namespace opensn
{

OpenSnRegisterObjectInNamespace(lbs, PointSource);

InputParameters
PointSource::GetInputParameters()
{
  InputParameters params;

  params.SetGeneralDescription("A multi-group isotropic point source.");
  params.SetClassName("Point Source");

  params.AddRequiredParameterArray("location", "The (x, y, z) coordinate of the point source.");
  params.AddOptionalParameterArray(
    "strength", std::vector<double>(), "The group-wise point source strength");
  params.AddOptionalParameter<std::shared_ptr<GroupTimeFunction>>(
    "strength_function",
    std::shared_ptr<GroupTimeFunction>{},
    "Function defining group-wise strengths as a function of (group, time).");
  params.AddOptionalParameter("start_time",
                              -std::numeric_limits<double>::infinity(),
                              "Time at which the source becomes active.");
  params.AddOptionalParameter("end_time",
                              std::numeric_limits<double>::infinity(),
                              "Time at which the source becomes inactive.");

  return params;
}

std::shared_ptr<PointSource>
PointSource::Create(const ParameterBlock& params)
{
  auto& factory = opensn::ObjectFactory::GetInstance();
  return factory.Create<PointSource>("lbs::PointSource", params);
}

PointSource::PointSource(const InputParameters& params)
  : location_(params.GetParamVectorValue<double>("location")),
    strength_(params.GetParamVectorValue<double>("strength")),
    strength_function_(params.GetSharedPtrParam<GroupTimeFunction>("strength_function", false)),
    start_time_(params.GetParamValue<double>("start_time")),
    end_time_(params.GetParamValue<double>("end_time"))
{
  const bool has_strength = not strength_.empty();
  const bool has_strength_func = static_cast<bool>(strength_function_);
  if ((has_strength ? 1 : 0) + (has_strength_func ? 1 : 0) == 0)
    throw std::invalid_argument("Either a strength vector or strength_function must be provided.");
  if ((has_strength ? 1 : 0) + (has_strength_func ? 1 : 0) > 1)
    throw std::invalid_argument(
      "Specify only one of strength or strength_function for a point source.");
  if (has_strength_func && (start_time_ != -std::numeric_limits<double>::infinity() ||
                            end_time_ != std::numeric_limits<double>::infinity()))
    throw std::invalid_argument("strength_function cannot be used with start_time/end_time. "
                                "Define time dependence in the callback or omit the time bounds.");
  if (not strength_.empty() &&
      std::all_of(strength_.begin(), strength_.end(), [](double x) { return x == 0.0; }))
    log.Log0Warning() << "Point source at " << location_.PrintStr() << " "
                      << "does not have a non-zero source strength.";
}

void
PointSource::Initialize(const LBSProblem& lbs_problem)
{
  if (not strength_function_)
  {
    OpenSnLogicalErrorIf(strength_.size() != lbs_problem.GetNumGroups(),
                         "Incompatible point source strength vector at location " +
                           location_.PrintStr() + ". " + "There are " +
                           std::to_string(lbs_problem.GetNumGroups()) + " energy groups, but " +
                           std::to_string(strength_.size()) + " source strength values.");
  }

  // Get info from solver
  const auto& grid = lbs_problem.GetGrid();
  const auto& discretization = lbs_problem.GetSpatialDiscretization();
  const auto& unit_cell_matrices = lbs_problem.GetUnitCellMatrices();
  const auto& ghost_unit_cell_matrices = lbs_problem.GetUnitGhostCellMatrices();

  // Find local subscribers
  double total_volume = 0.0;
  std::vector<Subscriber> subscribers;
  for (const auto& cell : grid->local_cells)
  {
    if (grid->CheckPointInsideCell(cell, location_))
    {
      const auto& cell_mapping = discretization.GetCellMapping(cell);
      const auto& fe_values = unit_cell_matrices[cell.local_id];

      // Map the point source to the finite element space
      Vector<double> shape_vals;
      cell_mapping.ShapeValues(location_, shape_vals);
      const auto M_inv = Inverse(fe_values.intV_shapeI_shapeJ);
      const auto node_wgts = Mult(M_inv, shape_vals);

      // Increment the total volume
      total_volume += cell.volume;

      // Add to subscribers
      subscribers.push_back(Subscriber{cell.volume, cell.local_id, shape_vals, node_wgts});
    }
  }

  // If the point source lies on a partition boundary, ghost cells must be
  // added to the total volume.
  auto ghost_global_ids = grid->cells.GetGhostGlobalIDs();
  for (uint64_t global_id : ghost_global_ids)
  {
    const auto& nbr_cell = grid->cells[global_id];
    if (grid->CheckPointInsideCell(nbr_cell, location_))
    {
      const auto& fe_values = ghost_unit_cell_matrices.at(nbr_cell.global_id);
      total_volume +=
        std::accumulate(fe_values.intV_shapeI.begin(), fe_values.intV_shapeI.end(), 0.0);
    }
  }

  // Create the actual subscriber list
  subscribers_.clear();
  for (const auto& sub : subscribers)
  {
    subscribers_.push_back(
      {sub.volume_weight / total_volume, sub.cell_local_id, sub.shape_values, sub.node_weights});

    std::stringstream ss;
    ss << "Point source at " << location_.PrintStr() << " assigned to cell "
       << grid->local_cells[sub.cell_local_id].global_id << " with shape values [ ";
    for (const auto& value : sub.shape_values)
      ss << value << " ";
    ss << "] and volume weight " << sub.volume_weight / total_volume;
    log.LogAll() << ss.str();
  }

  size_t num_local_subs = subscribers_.size();
  mpi_comm.all_reduce(num_local_subs, num_global_subscribers_, mpi::op::sum<size_t>());

  log.LogAll() << "Point source has " << num_local_subs << " subscribing cells on processor "
               << mpi_comm.rank();
  log.Log() << "Point source has " << num_global_subscribers_ << " global subscribing cells.";
}

bool
PointSource::IsActive(double time) const
{
  return time >= start_time_ && time <= end_time_;
}

std::vector<double>
PointSource::GetStrength(const double time, const unsigned int num_groups) const
{
  if (strength_function_)
  {
    std::vector<double> values(num_groups, 0.0);
    for (unsigned int g = 0; g < num_groups; ++g)
      values[g] = (*strength_function_)(g, time);
    return values;
  }
  return strength_;
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "modules/linear_boltzmann_solvers/lbs_problem/point_source/point_source.h"
5	#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
6	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
7	#include "framework/math/functions/function.h"
8	#include "framework/object_factory.h"
9	#include "framework/logging/log.h"
10	#include "framework/runtime.h"
11	#include <numeric>
12	#include <limits>
13
14	namespace opensn
15	{
16
17	OpenSnRegisterObjectInNamespace(lbs, PointSource);
18
19	InputParameters
20	PointSource::GetInputParameters()	9✔
21	{
22	InputParameters params;	9✔
23
24	params.SetGeneralDescription("A multi-group isotropic point source.");	18✔
25	params.SetClassName("Point Source");	18✔
26
27	params.AddRequiredParameterArray("location", "The (x, y, z) coordinate of the point source.");	18✔
28	params.AddOptionalParameterArray(	18✔
29	"strength", std::vector<double>(), "The group-wise point source strength");	18✔
30	params.AddOptionalParameter<std::shared_ptr<GroupTimeFunction>>(	18✔
31	"strength_function",
32	std::shared_ptr<GroupTimeFunction>{},	18✔
33	"Function defining group-wise strengths as a function of (group, time).");
34	params.AddOptionalParameter("start_time",	18✔
35	-std::numeric_limits<double>::infinity(),
36	"Time at which the source becomes active.");
37	params.AddOptionalParameter("end_time",	18✔
38	std::numeric_limits<double>::infinity(),
39	"Time at which the source becomes inactive.");
40
41	return params;	9✔
UNCOV 42	}	×
43
44	std::shared_ptr<PointSource>
45	PointSource::Create(const ParameterBlock& params)	9✔
46	{
47	auto& factory = opensn::ObjectFactory::GetInstance();	9✔
48	return factory.Create<PointSource>("lbs::PointSource", params);	18✔
49	}
50
51	PointSource::PointSource(const InputParameters& params)	9✔
52	: location_(params.GetParamVectorValue<double>("location")),	9✔
53	strength_(params.GetParamVectorValue<double>("strength")),	9✔
54	strength_function_(params.GetSharedPtrParam<GroupTimeFunction>("strength_function", false)),	9✔
55	start_time_(params.GetParamValue<double>("start_time")),	9✔
56	end_time_(params.GetParamValue<double>("end_time"))	9✔
57	{
58	const bool has_strength = not strength_.empty();	9✔
59	const bool has_strength_func = static_cast<bool>(strength_function_);	9✔
60	if ((has_strength ? 1 : 0) + (has_strength_func ? 1 : 0) == 0)	9✔
UNCOV 61	throw std::invalid_argument("Either a strength vector or strength_function must be provided.");	×
62	if ((has_strength ? 1 : 0) + (has_strength_func ? 1 : 0) > 1)	9✔
UNCOV 63	throw std::invalid_argument(	×
UNCOV 64	"Specify only one of strength or strength_function for a point source.");	×
65	if (has_strength_func && (start_time_ != -std::numeric_limits<double>::infinity() \|\|	9✔
UNCOV 66	end_time_ != std::numeric_limits<double>::infinity()))	×
UNCOV 67	throw std::invalid_argument("strength_function cannot be used with start_time/end_time. "	×
UNCOV 68	"Define time dependence in the callback or omit the time bounds.");	×
69	if (not strength_.empty() &&	9✔
70	std::all_of(strength_.begin(), strength_.end(), [](double x) { return x == 0.0; }))	9✔
UNCOV 71	log.Log0Warning() << "Point source at " << location_.PrintStr() << " "	×
UNCOV 72	<< "does not have a non-zero source strength.";	×
73	}	9✔
74
75	void
76	PointSource::Initialize(const LBSProblem& lbs_problem)	17✔
77	{
78	if (not strength_function_)	17✔
79	{
80	OpenSnLogicalErrorIf(strength_.size() != lbs_problem.GetNumGroups(),	17✔
81	"Incompatible point source strength vector at location " +
82	location_.PrintStr() + ". " + "There are " +
83	std::to_string(lbs_problem.GetNumGroups()) + " energy groups, but " +
84	std::to_string(strength_.size()) + " source strength values.");
85	}
86
87	// Get info from solver
88	const auto& grid = lbs_problem.GetGrid();	17✔
89	const auto& discretization = lbs_problem.GetSpatialDiscretization();	17✔
90	const auto& unit_cell_matrices = lbs_problem.GetUnitCellMatrices();	17✔
91	const auto& ghost_unit_cell_matrices = lbs_problem.GetUnitGhostCellMatrices();	17✔
92
93	// Find local subscribers
94	double total_volume = 0.0;	17✔
95	std::vector<Subscriber> subscribers;	17✔
96	for (const auto& cell : grid->local_cells)	14,538✔
97	{
98	if (grid->CheckPointInsideCell(cell, location_))	14,521✔
99	{
100	const auto& cell_mapping = discretization.GetCellMapping(cell);	5✔
101	const auto& fe_values = unit_cell_matrices[cell.local_id];	5✔
102
103	// Map the point source to the finite element space
104	Vector<double> shape_vals;	5✔
105	cell_mapping.ShapeValues(location_, shape_vals);	5✔
106	const auto M_inv = Inverse(fe_values.intV_shapeI_shapeJ);	5✔
107	const auto node_wgts = Mult(M_inv, shape_vals);	5✔
108
109	// Increment the total volume
110	total_volume += cell.volume;	5✔
111
112	// Add to subscribers
113	subscribers.push_back(Subscriber{cell.volume, cell.local_id, shape_vals, node_wgts});	5✔
114	}	15✔
115	}
116
117	// If the point source lies on a partition boundary, ghost cells must be
118	// added to the total volume.
119	auto ghost_global_ids = grid->cells.GetGhostGlobalIDs();	17✔
120	for (uint64_t global_id : ghost_global_ids)	1,025✔
121	{
122	const auto& nbr_cell = grid->cells[global_id];	1,008✔
123	if (grid->CheckPointInsideCell(nbr_cell, location_))	1,008✔
124	{
UNCOV 125	const auto& fe_values = ghost_unit_cell_matrices.at(nbr_cell.global_id);	×
UNCOV 126	total_volume +=	×
UNCOV 127	std::accumulate(fe_values.intV_shapeI.begin(), fe_values.intV_shapeI.end(), 0.0);	×
128	}
129	}
130
131	// Create the actual subscriber list
132	subscribers_.clear();	17✔
133	for (const auto& sub : subscribers)	22✔
134	{
135	subscribers_.push_back(	5✔
136	{sub.volume_weight / total_volume, sub.cell_local_id, sub.shape_values, sub.node_weights});	10✔
137
138	std::stringstream ss;	5✔
139	ss << "Point source at " << location_.PrintStr() << " assigned to cell "	10✔
140	<< grid->local_cells[sub.cell_local_id].global_id << " with shape values [ ";	5✔
141	for (const auto& value : sub.shape_values)	25✔
142	ss << value << " ";	20✔
143	ss << "] and volume weight " << sub.volume_weight / total_volume;	5✔
144	log.LogAll() << ss.str();	10✔
145	}	5✔
146
147	size_t num_local_subs = subscribers_.size();	17✔
148	mpi_comm.all_reduce(num_local_subs, num_global_subscribers_, mpi::op::sum<size_t>());	17✔
149
150	log.LogAll() << "Point source has " << num_local_subs << " subscribing cells on processor "	17✔
151	<< mpi_comm.rank();	34✔
152	log.Log() << "Point source has " << num_global_subscribers_ << " global subscribing cells.";	51✔
153	}	39✔
154
155	bool
156	PointSource::IsActive(double time) const	17✔
157	{
158	return time >= start_time_ && time <= end_time_;	17✔
159	}
160
161	std::vector<double>
162	PointSource::GetStrength(const double time, const unsigned int num_groups) const	7✔
163	{
164	if (strength_function_)	7✔
165	{
UNCOV 166	std::vector<double> values(num_groups, 0.0);	×
UNCOV 167	for (unsigned int g = 0; g < num_groups; ++g)	×
UNCOV 168	values[g] = (*strength_function_)(g, time);	×
169	return values;
UNCOV 170	}	×
171	return strength_;	7✔
172	}
173
174	} // namespace opensn

Open-Sn / opensn / 21892714514

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