25903849510

Committed 14 May 2026 09:05PM UTC coverage: 76.114% (+0.4%) from 75.683%

Build # 25903849510

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1071 from wdhawkins/input_wildcard

Add wildcard support to regression test script.

Coverage Stats

22612 of 29708 relevant lines covered (76.11%)

65089871.83 hits per line

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

96.21

/modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.cc

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

#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/source_functions/source_function.h"
#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
#include "framework/mesh/mesh_continuum/mesh_continuum.h"
#include "framework/runtime.h"
#include "framework/logging/log.h"

namespace opensn
{

SourceFunction::SourceFunction(const LBSProblem& lbs_problem) : lbs_problem_(lbs_problem)
{
}

void
SourceFunction::operator()(const LBSGroupset& groupset,
                           std::vector<double>& q,
                           const std::vector<double>& phi,
                           const SourceFlags source_flags)
{
  if (source_flags.Empty())
    return;

  apply_fixed_src_ = (source_flags & APPLY_FIXED_SOURCES);
  apply_wgs_scatter_src_ = (source_flags & APPLY_WGS_SCATTER_SOURCES);
  apply_ags_scatter_src_ = (source_flags & APPLY_AGS_SCATTER_SOURCES);
  apply_wgs_fission_src_ = (source_flags & APPLY_WGS_FISSION_SOURCES);
  apply_ags_fission_src_ = (source_flags & APPLY_AGS_FISSION_SOURCES);
  suppress_wg_scatter_src_ = (source_flags & SUPPRESS_WG_SCATTER);

  // Get group setup
  gs_i_ = groupset.first_group;
  gs_f_ = groupset.last_group;

  first_grp_ = 0;
  last_grp_ = lbs_problem_.GetNumGroups() - 1;

  default_zero_src_.assign(lbs_problem_.GetNumGroups(), 0.0);

  const auto& cell_transport_views = lbs_problem_.GetCellTransportViews();

  const auto num_moments = lbs_problem_.GetNumMoments();
  const auto& ext_src_moments_local = lbs_problem_.GetExtSrcMomentsLocal();

  const auto& m_to_ell_em_map = groupset.quadrature->GetMomentToHarmonicsIndexMap();

  // Apply all nodal sources
  const auto& grid = lbs_problem_.GetGrid();
  for (const auto& cell : grid->local_cells)
  {
    const auto& transport_view = cell_transport_views[cell.local_id];
    cell_volume_ = transport_view.GetVolume();

    // Obtain xs
    const auto& xs = transport_view.GetXS();

    const auto& S = xs.GetTransferMatrices();
    const auto& F = xs.GetProductionMatrix();
    const auto& precursors = xs.GetPrecursors();
    const auto& nu_delayed_sigma_f = xs.GetNuDelayedSigmaF();

    // Loop over nodes
    const auto num_nodes = transport_view.GetNumNodes();
    for (int i = 0; i < num_nodes; ++i)
    {
      // Loop over moments
      for (unsigned int m = 0; m < num_moments; ++m)
      {
        const auto ell = m_to_ell_em_map[m].ell;
        const auto uk_map = transport_view.MapDOF(i, m, 0);
        const double* phi_im = &phi[uk_map];

        // Declare moment src
        fixed_src_moments_ = default_zero_src_.data();

        if (lbs_problem_.GetOptions().use_src_moments)
          fixed_src_moments_ = &ext_src_moments_local[uk_map];

        // Loop over groupset groups
        for (auto g = gs_i_; g <= gs_f_; ++g)
        {
          g_ = g;

          double rhs = 0.0;

          // Apply fixed sources
          if (apply_fixed_src_)
            rhs += FixedSourceMoments();

          // Apply scattering sources
          if (ell < S.size())
          {
            const auto& S_ell = S[ell];
            // Add Across GroupSet Scattering (AGS)
            if (apply_ags_scatter_src_)
              for (const auto& [_, gp, sigma_sm] : S_ell.Row(g))
                if (gp < gs_i_ or gp > gs_f_)
                  rhs += sigma_sm * phi_im[gp];

            // Add Within GroupSet Scattering (WGS)
            if (apply_wgs_scatter_src_)
              for (const auto& [_, gp, sigma_sm] : S_ell.Row(g))
                if (gp >= gs_i_ and gp <= gs_f_)
                {
                  if (suppress_wg_scatter_src_ and g_ == gp)
                    continue;
                  rhs += sigma_sm * phi_im[gp];
                }
          }

          // Apply fission sources
          if (xs.IsFissionable() and ell == 0)
          {
            const auto& F_g = F[g];
            if (apply_ags_fission_src_)
              for (size_t gp = first_grp_; gp <= last_grp_; ++gp)
                if (gp < gs_i_ or gp > gs_f_)
                  rhs += F_g[gp] * phi_im[gp];

            if (apply_wgs_fission_src_)
              for (size_t gp = gs_i_; gp <= gs_f_; ++gp)
                rhs += F_g[gp] * phi_im[gp];

            if (lbs_problem_.GetOptions().use_precursors)
              rhs += DelayedFission(precursors, nu_delayed_sigma_f, &phi[uk_map]);
          }

          // Add to destination vector
          q[uk_map + g] += rhs;

        } // for g
      } // for m
    } // for dof i
  } // for cell

  AddAdditionalSources(groupset, q, phi, source_flags);
}

double
SourceFunction::FixedSourceMoments() const
{
  return fixed_src_moments_[g_];
}

double
SourceFunction::DelayedFission(const PrecursorList& precursors,
                               const std::vector<double>& nu_delayed_sigma_f,
                               const double* phi) const
{
  double value = 0.0;
  if (apply_ags_fission_src_)
    for (auto gp = first_grp_; gp <= last_grp_; ++gp)
      if (gp < gs_i_ or gp > gs_f_)
        for (const auto& precursor : precursors)
          value += precursor.emission_spectrum[g_] * precursor.fractional_yield *
                   nu_delayed_sigma_f[gp] * phi[gp];

  if (apply_wgs_fission_src_)
    for (auto gp = gs_i_; gp <= gs_f_; ++gp)
      for (const auto& precursor : precursors)
        value += precursor.emission_spectrum[g_] * precursor.fractional_yield *
                 nu_delayed_sigma_f[gp] * phi[gp];

  return value;
}

void
SourceFunction::AddPointSources(const LBSGroupset& groupset,
                                std::vector<double>& q,
                                const std::vector<double>& /*unused*/,
                                const SourceFlags source_flags)
{
  const bool apply_fixed_src = (source_flags & APPLY_FIXED_SOURCES);

  const auto& transport_views = lbs_problem_.GetCellTransportViews();

  const auto gs_i = groupset.first_group;
  const auto gs_f = groupset.last_group;
  const double source_time = lbs_problem_.GetTime();

  // Apply point sources
  if (not lbs_problem_.GetOptions().use_src_moments and apply_fixed_src)
  {
    for (const auto& point_source : lbs_problem_.GetPointSources())
    {
      if (not point_source->IsActive(source_time))
        continue;

      for (const auto& subscriber : point_source->GetSubscribers())
      {
        const auto& transport_view = transport_views[subscriber.cell_local_id];

        const auto strength = point_source->GetStrength(source_time, lbs_problem_.GetNumGroups());
        const auto& node_weights = subscriber.node_weights;
        const auto volume_weight = subscriber.volume_weight;

        for (int i = 0; i < transport_view.GetNumNodes(); ++i)
        {
          const auto uk_map = transport_view.MapDOF(i, 0, 0);
          for (auto g = gs_i; g <= gs_f; ++g)
            q[uk_map + g] += strength[g] * node_weights(i) * volume_weight;
        } // for node i
      } // for subscriber
    } // for point source
  }
}

void
SourceFunction::AddVolumetricSources(const LBSGroupset& groupset,
                                     std::vector<double>& q,
                                     const std::vector<double>& phi,
                                     const SourceFlags source_flags)
{
  const bool apply_fixed_src = source_flags & APPLY_FIXED_SOURCES;

  const auto& grid = lbs_problem_.GetGrid();
  const auto& discretization = lbs_problem_.GetSpatialDiscretization();
  const auto& cell_transport_views = lbs_problem_.GetCellTransportViews();
  const auto num_groups = lbs_problem_.GetNumGroups();

  const auto gs_i = groupset.first_group;
  const auto gs_f = groupset.last_group;
  const double source_time = lbs_problem_.GetTime();

  // Go through each volumetric source, and its subscribing cells
  if (not lbs_problem_.GetOptions().use_src_moments and apply_fixed_src)
  {
    for (const auto& volumetric_source : lbs_problem_.GetVolumetricSources())
    {
      if (not volumetric_source->IsActive(source_time))
        continue;

      for (const auto local_id : volumetric_source->GetSubscribers())
      {
        const auto& cell = grid->local_cells[local_id];
        const auto& transport_view = cell_transport_views[local_id];
        const auto nodes = discretization.GetCellNodeLocations(cell);
        const auto num_cell_nodes = discretization.GetCellNumNodes(cell);

        // Go through each of the cell nodes
        for (size_t i = 0; i < num_cell_nodes; ++i)
        {
          // Compute group-wise values for this node
          const auto src = volumetric_source->Evaluate(cell, nodes[i], num_groups, source_time);

          // Contribute to the source moments
          const auto dof_map = transport_view.MapDOF(i, 0, 0);
          for (auto g = gs_i; g <= gs_f; ++g)
            q[dof_map + g] += src[g];
        } // for node i
      } // for subscriber
    } // for volumetric source
  }
}

} // 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/discrete_ordinates_problem/source_functions/source_function.h"
5	#include "modules/linear_boltzmann_solvers/lbs_problem/lbs_problem.h"
6	#include "framework/mesh/mesh_continuum/mesh_continuum.h"
7	#include "framework/runtime.h"
8	#include "framework/logging/log.h"
9
10	namespace opensn
11	{
12
13	SourceFunction::SourceFunction(const LBSProblem& lbs_problem) : lbs_problem_(lbs_problem)	1,196✔
14	{
15	}	1,196✔
16
17	void
18	SourceFunction::operator()(const LBSGroupset& groupset,	114,803✔
19	std::vector<double>& q,
20	const std::vector<double>& phi,
21	const SourceFlags source_flags)
22	{
23	if (source_flags.Empty())	114,803✔
24	return;	×
25
26	apply_fixed_src_ = (source_flags & APPLY_FIXED_SOURCES);	114,803✔
27	apply_wgs_scatter_src_ = (source_flags & APPLY_WGS_SCATTER_SOURCES);	114,803✔
28	apply_ags_scatter_src_ = (source_flags & APPLY_AGS_SCATTER_SOURCES);	114,803✔
29	apply_wgs_fission_src_ = (source_flags & APPLY_WGS_FISSION_SOURCES);	114,803✔
30	apply_ags_fission_src_ = (source_flags & APPLY_AGS_FISSION_SOURCES);	114,803✔
31	suppress_wg_scatter_src_ = (source_flags & SUPPRESS_WG_SCATTER);	114,803✔
32
33	// Get group setup
34	gs_i_ = groupset.first_group;	114,803✔
35	gs_f_ = groupset.last_group;	114,803✔
36
37	first_grp_ = 0;	114,803✔
38	last_grp_ = lbs_problem_.GetNumGroups() - 1;	114,803✔
39
40	default_zero_src_.assign(lbs_problem_.GetNumGroups(), 0.0);	114,803✔
41
42	const auto& cell_transport_views = lbs_problem_.GetCellTransportViews();	114,803✔
43
44	const auto num_moments = lbs_problem_.GetNumMoments();	114,803✔
45	const auto& ext_src_moments_local = lbs_problem_.GetExtSrcMomentsLocal();	114,803✔
46
47	const auto& m_to_ell_em_map = groupset.quadrature->GetMomentToHarmonicsIndexMap();	114,803✔
48
49	// Apply all nodal sources
50	const auto& grid = lbs_problem_.GetGrid();	114,803✔
51	for (const auto& cell : grid->local_cells)	59,560,638✔
52	{
53	const auto& transport_view = cell_transport_views[cell.local_id];	59,445,835✔
54	cell_volume_ = transport_view.GetVolume();	59,445,835✔
55
56	// Obtain xs
57	const auto& xs = transport_view.GetXS();	59,445,835✔
58
59	const auto& S = xs.GetTransferMatrices();	59,445,835✔
60	const auto& F = xs.GetProductionMatrix();	59,445,835✔
61	const auto& precursors = xs.GetPrecursors();	59,445,835✔
62	const auto& nu_delayed_sigma_f = xs.GetNuDelayedSigmaF();	59,445,835✔
63
64	// Loop over nodes
65	const auto num_nodes = transport_view.GetNumNodes();	59,445,835✔
66	for (int i = 0; i < num_nodes; ++i)	320,058,165✔
67	{
68	// Loop over moments
69	for (unsigned int m = 0; m < num_moments; ++m)	1,167,411,736✔
70	{
71	const auto ell = m_to_ell_em_map[m].ell;	906,799,406✔
72	const auto uk_map = transport_view.MapDOF(i, m, 0);	906,799,406✔
73	const double* phi_im = &phi[uk_map];	906,799,406✔
74
75	// Declare moment src
76	fixed_src_moments_ = default_zero_src_.data();	906,799,406✔
77
78	if (lbs_problem_.GetOptions().use_src_moments)	906,799,406✔
79	fixed_src_moments_ = &ext_src_moments_local[uk_map];	1,572,864✔
80
81	// Loop over groupset groups
82	for (auto g = gs_i_; g <= gs_f_; ++g)	2,147,483,647✔
83	{
84	g_ = g;	2,147,483,647✔
85
86	double rhs = 0.0;	2,147,483,647✔
87
88	// Apply fixed sources
89	if (apply_fixed_src_)	2,147,483,647✔
90	rhs += FixedSourceMoments();	2,147,483,647✔
91
92	// Apply scattering sources
93	if (ell < S.size())	2,147,483,647✔
94	{
95	const auto& S_ell = S[ell];	2,147,483,647✔
96	// Add Across GroupSet Scattering (AGS)
97	if (apply_ags_scatter_src_)	2,147,483,647✔
98	for (const auto& [_, gp, sigma_sm] : S_ell.Row(g))	2,147,483,647✔
99	if (gp < gs_i_ or gp > gs_f_)	2,147,483,647✔
100	rhs += sigma_sm * phi_im[gp];	301,364,320✔
101
102	// Add Within GroupSet Scattering (WGS)
103	if (apply_wgs_scatter_src_)	2,147,483,647✔
104	for (const auto& [_, gp, sigma_sm] : S_ell.Row(g))	2,147,483,647✔
105	if (gp >= gs_i_ and gp <= gs_f_)	2,147,483,647✔
106	{
107	if (suppress_wg_scatter_src_ and g_ == gp)	2,147,483,647✔
108	continue;	36,020,056✔
109	rhs += sigma_sm * phi_im[gp];	2,147,483,647✔
110	}
111	}
112
113	// Apply fission sources
114	if (xs.IsFissionable() and ell == 0)	2,147,483,647✔
115	{
116	const auto& F_g = F[g];	436,877,208✔
117	if (apply_ags_fission_src_)	436,877,208✔
118	for (size_t gp = first_grp_; gp <= last_grp_; ++gp)	2,147,483,647✔
119	if (gp < gs_i_ or gp > gs_f_)	2,147,483,647✔
120	rhs += F_g[gp] * phi_im[gp];	12,800✔
121
122	if (apply_wgs_fission_src_)	436,877,208✔
123	for (size_t gp = gs_i_; gp <= gs_f_; ++gp)	2,147,483,647✔
124	rhs += F_g[gp] * phi_im[gp];	2,147,483,647✔
125
126	if (lbs_problem_.GetOptions().use_precursors)	436,877,208✔
127	rhs += DelayedFission(precursors, nu_delayed_sigma_f, &phi[uk_map]);	345,528,248✔
128	}
129
130	// Add to destination vector
131	q[uk_map + g] += rhs;	2,147,483,647✔
132
133	} // for g
134	} // for m
135	} // for dof i
136	} // for cell
137
138	AddAdditionalSources(groupset, q, phi, source_flags);	114,803✔
139	}	114,803✔
140
141	double
142	SourceFunction::FixedSourceMoments() const	2,147,483,647✔
143	{
144	return fixed_src_moments_[g_];	2,147,483,647✔
145	}
146
147	double
148	SourceFunction::DelayedFission(const PrecursorList& precursors,	345,173,000✔
149	const std::vector<double>& nu_delayed_sigma_f,
150	const double* phi) const
151	{
152	double value = 0.0;	345,173,000✔
153	if (apply_ags_fission_src_)	345,173,000✔
154	for (auto gp = first_grp_; gp <= last_grp_; ++gp)	2,147,483,647✔
155	if (gp < gs_i_ or gp > gs_f_)	2,147,483,647✔
156	for (const auto& precursor : precursors)	×
157	value += precursor.emission_spectrum[g_] * precursor.fractional_yield *	×
158	nu_delayed_sigma_f[gp] * phi[gp];	×
159
160	if (apply_wgs_fission_src_)	345,173,000✔
161	for (auto gp = gs_i_; gp <= gs_f_; ++gp)	2,147,483,647✔
162	for (const auto& precursor : precursors)	2,147,483,647✔
163	value += precursor.emission_spectrum[g_] * precursor.fractional_yield *	306,308✔
164	nu_delayed_sigma_f[gp] * phi[gp];	306,308✔
165
166	return value;	345,173,000✔
167	}
168
169	void
170	SourceFunction::AddPointSources(const LBSGroupset& groupset,	114,803✔
171	std::vector<double>& q,
172	const std::vector<double>& /unused/,
173	const SourceFlags source_flags)
174	{
175	const bool apply_fixed_src = (source_flags & APPLY_FIXED_SOURCES);	114,803✔
176
177	const auto& transport_views = lbs_problem_.GetCellTransportViews();	114,803✔
178
179	const auto gs_i = groupset.first_group;	114,803✔
180	const auto gs_f = groupset.last_group;	114,803✔
181	const double source_time = lbs_problem_.GetTime();	114,803✔
182
183	// Apply point sources
184	if (not lbs_problem_.GetOptions().use_src_moments and apply_fixed_src)	114,803✔
185	{
186	for (const auto& point_source : lbs_problem_.GetPointSources())	27,650✔
187	{
188	if (not point_source->IsActive(source_time))	16✔
189	continue;	×
190
191	for (const auto& subscriber : point_source->GetSubscribers())	20✔
192	{
193	const auto& transport_view = transport_views[subscriber.cell_local_id];	4✔
194
195	const auto strength = point_source->GetStrength(source_time, lbs_problem_.GetNumGroups());	4✔
196	const auto& node_weights = subscriber.node_weights;	4✔
197	const auto volume_weight = subscriber.volume_weight;	4✔
198
199	for (int i = 0; i < transport_view.GetNumNodes(); ++i)	20✔
200	{
201	const auto uk_map = transport_view.MapDOF(i, 0, 0);	16✔
202	for (auto g = gs_i; g <= gs_f; ++g)	104✔
203	q[uk_map + g] += strength[g] * node_weights(i) * volume_weight;	88✔
204	} // for node i
205	} // for subscriber	4✔
206	} // for point source
207	}
208	}	114,803✔
209
210	void
211	SourceFunction::AddVolumetricSources(const LBSGroupset& groupset,	114,803✔
212	std::vector<double>& q,
213	const std::vector<double>& phi,
214	const SourceFlags source_flags)
215	{
216	const bool apply_fixed_src = source_flags & APPLY_FIXED_SOURCES;	114,803✔
217
218	const auto& grid = lbs_problem_.GetGrid();	114,803✔
219	const auto& discretization = lbs_problem_.GetSpatialDiscretization();	114,803✔
220	const auto& cell_transport_views = lbs_problem_.GetCellTransportViews();	114,803✔
221	const auto num_groups = lbs_problem_.GetNumGroups();	114,803✔
222
223	const auto gs_i = groupset.first_group;	114,803✔
224	const auto gs_f = groupset.last_group;	114,803✔
225	const double source_time = lbs_problem_.GetTime();	114,803✔
226
227	// Go through each volumetric source, and its subscribing cells
228	if (not lbs_problem_.GetOptions().use_src_moments and apply_fixed_src)	114,803✔
229	{
230	for (const auto& volumetric_source : lbs_problem_.GetVolumetricSources())	35,710✔
231	{
232	if (not volumetric_source->IsActive(source_time))	8,076✔
233	continue;	472✔
234
235	for (const auto local_id : volumetric_source->GetSubscribers())	5,308,170✔
236	{
237	const auto& cell = grid->local_cells[local_id];	5,300,566✔
238	const auto& transport_view = cell_transport_views[local_id];	5,300,566✔
239	const auto nodes = discretization.GetCellNodeLocations(cell);	5,300,566✔
240	const auto num_cell_nodes = discretization.GetCellNumNodes(cell);	5,300,566✔
241
242	// Go through each of the cell nodes
243	for (size_t i = 0; i < num_cell_nodes; ++i)	33,659,858✔
244	{
245	// Compute group-wise values for this node
246	const auto src = volumetric_source->Evaluate(cell, nodes[i], num_groups, source_time);	28,359,292✔
247
248	// Contribute to the source moments
249	const auto dof_map = transport_view.MapDOF(i, 0, 0);	28,359,292✔
250	for (auto g = gs_i; g <= gs_f; ++g)	141,387,860✔
251	q[dof_map + g] += src[g];	113,028,568✔
252	} // for node i	28,359,292✔
253	} // for subscriber	5,300,566✔
254	} // for volumetric source
255	}
256	}	114,803✔
257
258	} // namespace opensn

Open-Sn / opensn / 25903849510

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