20201962283

Committed 12 Dec 2025 06:55PM UTC coverage: 74.333%. Remained the same

Build # 20201962283

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #859 from wdhawkins/td_source_driver

Adding time-dependent solver and time-dependent sources

Run Details

367 of 398 new or added lines in 23 files covered. (92.21%)

359 existing lines in 17 files now uncovered.

18610 of 25036 relevant lines covered (74.33%)

68925445.1 hits per line

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

95.52

/modules/linear_boltzmann_solvers/lbs_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/lbs_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"
#include "caliper/cali.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)
{
  CALI_CXX_MARK_SCOPE("SourceFunction::operator");

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

  const auto& densities = lbs_problem_.GetDensitiesLocal();

  // Get group setup
  gs_i_ = static_cast<size_t>(groupset.groups.front().id);
  gs_f_ = static_cast<size_t>(groupset.groups.back().id);

  first_grp_ = static_cast<size_t>(lbs_problem_.GetGroups().front().id);
  last_grp_ = static_cast<size_t>(lbs_problem_.GetGroups().back().id);

  default_zero_src_.assign(lbs_problem_.GetGroups().size(), 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& rho = densities[cell.local_id];
    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 (size_t 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 (size_t 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 += rho * 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 += rho * 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 += rho * F_g[gp] * phi_im[gp];

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

            if (lbs_problem_.GetOptions().use_precursors)
              rhs += DelayedFission(precursors, rho, 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 double& rho,
                               const std::vector<double>& nu_delayed_sigma_f,
                               const double* phi) const
{
  double value = 0.0;
  if (apply_ags_fission_src_)
    for (size_t 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 * rho *
                   nu_delayed_sigma_f[gp] * phi[gp];

  if (apply_wgs_fission_src_)
    for (size_t gp = gs_i_; gp <= gs_f_; ++gp)
      for (const auto& precursor : precursors)
        value += precursor.emission_spectrum[g_] * precursor.fractional_yield * rho *
                 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.groups.front().id;
  const auto gs_f = groupset.groups.back().id;
  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();
        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 (int 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.groups.front().id;
  const auto gs_f = groupset.groups.back().id;
  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)(cell, nodes[i], num_groups);

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

Open-Sn / opensn / 20201962283

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