21887062786

Committed 10 Feb 2026 11:55PM UTC coverage: 81.899% (-0.1%) from 82.009%

Build # 21887062786

Build Type

Pull #3493

github

Committed by

web-flow

Commit Message

Merge 320b68711 into 3f20a5e22

Pull Request Pull Request #3493: Implement vector fitting to replace external `vectfit` package

Run Details

17361 of 24292 branches covered (71.47%)

Branch coverage included in aggregate %.

185 of 218 new or added lines in 3 files covered. (84.86%)

2770 existing lines in 69 files now uncovered.

56369 of 65733 relevant lines covered (85.75%)

51144917.16 hits per line

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

88.41

/src/random_ray/linear_source_domain.cpp

#include "openmc/random_ray/linear_source_domain.h"

#include "openmc/cell.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/output.h"
#include "openmc/plot.h"
#include "openmc/random_ray/random_ray.h"
#include "openmc/simulation.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/timer.h"

namespace openmc {

//==============================================================================
// LinearSourceDomain implementation
//==============================================================================

void LinearSourceDomain::batch_reset()
{
  FlatSourceDomain::batch_reset();
#pragma omp parallel for
  for (int64_t sr = 0; sr < n_source_regions(); sr++) {
    source_regions_.centroid_iteration(sr) = {0.0, 0.0, 0.0};
    source_regions_.mom_matrix(sr) = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  }
#pragma omp parallel for
  for (int64_t se = 0; se < n_source_elements(); se++) {
    source_regions_.flux_moments_new(se) = {0.0, 0.0, 0.0};
  }
}

void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)
{
  // Reset all source regions to zero (important for void regions)
  for (int g = 0; g < negroups_; g++) {
    srh.source(g) = 0.0;
  }

  // Add scattering + fission source
  int material = srh.material();
  int temp = srh.temperature_idx();
  double density_mult = srh.density_mult();
  if (material != MATERIAL_VOID) {
    double inverse_k_eff = 1.0 / k_eff_;
    MomentMatrix invM = srh.mom_matrix().inverse();

    for (int g_out = 0; g_out < negroups_; g_out++) {
      double sigma_t =
        sigma_t_[(material * ntemperature_ + temp) * negroups_ + g_out] *
        density_mult;

      double scatter_flat = 0.0f;
      double fission_flat = 0.0f;
      MomentArray scatter_linear = {0.0, 0.0, 0.0};
      MomentArray fission_linear = {0.0, 0.0, 0.0};

      for (int g_in = 0; g_in < negroups_; g_in++) {
        // Handles for the flat and linear components of the flux
        double flux_flat = srh.scalar_flux_old(g_in);
        MomentArray flux_linear = srh.flux_moments_old(g_in);

        // Handles for cross sections
        double sigma_s =
          sigma_s_[((material * ntemperature_ + temp) * negroups_ + g_out) *
                     negroups_ +
                   g_in] *
          density_mult;
        double nu_sigma_f =
          nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g_in] *
          density_mult;
        double chi =
          chi_[(material * ntemperature_ + temp) * negroups_ + g_out];

        // Compute source terms for flat and linear components of the flux
        scatter_flat += sigma_s * flux_flat;
        scatter_linear += sigma_s * flux_linear;
        if (settings::create_fission_neutrons) {
          fission_flat += nu_sigma_f * flux_flat * chi;
          fission_linear += nu_sigma_f * flux_linear * chi;
        }
      }

      // Compute the flat source term
      srh.source(g_out) =
        (scatter_flat + fission_flat * inverse_k_eff) / sigma_t;

      // Compute the linear source terms. In the first 10 iterations when the
      // centroids and spatial moments are not well known, we will leave the
      // source gradients as zero so as to avoid causing any numerical
      // instability. If a negative source is encountered, this region must be
      // very small/noisy or have poorly developed spatial moments, so we zero
      // the source gradients (effectively making this a flat source region
      // temporarily), so as to improve stability.
      if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) {
        srh.source_gradients(g_out) =
          invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
      } else {
        srh.source_gradients(g_out) = {0.0, 0.0, 0.0};
      }
    }
  }

  // Add external source if in fixed source mode
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
    for (int g = 0; g < negroups_; g++) {
      srh.source(g) += srh.external_source(g);
    }
  }
}

void LinearSourceDomain::normalize_scalar_flux_and_volumes(
  double total_active_distance_per_iteration)
{
  double normalization_factor = 1.0 / total_active_distance_per_iteration;
  double volume_normalization_factor =
    1.0 / (total_active_distance_per_iteration * simulation::current_batch);

// Normalize flux to total distance travelled by all rays this iteration
#pragma omp parallel for
  for (int64_t se = 0; se < n_source_elements(); se++) {
    source_regions_.scalar_flux_new(se) *= normalization_factor;
    source_regions_.flux_moments_new(se) *= normalization_factor;
  }

// Accumulate cell-wise ray length tallies collected this iteration, then
// update the simulation-averaged cell-wise volume estimates
#pragma omp parallel for
  for (int64_t sr = 0; sr < n_source_regions(); sr++) {
    source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr);
    source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr);
    source_regions_.volume_t(sr) += source_regions_.volume(sr);
    source_regions_.volume_sq_t(sr) += source_regions_.volume_sq(sr);
    source_regions_.volume_naive(sr) =
      source_regions_.volume(sr) * normalization_factor;
    source_regions_.volume(sr) =
      source_regions_.volume_t(sr) * volume_normalization_factor;
    source_regions_.volume_sq(sr) =
      source_regions_.volume_sq_t(sr) / source_regions_.volume_t(sr);
    if (source_regions_.volume_t(sr) > 0.0) {
      double inv_volume = 1.0 / source_regions_.volume_t(sr);
      source_regions_.centroid(sr) = source_regions_.centroid_t(sr);
      source_regions_.centroid(sr) *= inv_volume;
      source_regions_.mom_matrix(sr) = source_regions_.mom_matrix_t(sr);
      source_regions_.mom_matrix(sr) *= inv_volume;
    }
  }
}

void LinearSourceDomain::set_flux_to_flux_plus_source(
  int64_t sr, double volume, int g)
{
  int material = source_regions_.material(sr);
  if (material == MATERIAL_VOID) {
    FlatSourceDomain::set_flux_to_flux_plus_source(sr, volume, g);
  } else {
    source_regions_.scalar_flux_new(sr, g) /= volume;
    source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);
  }
  // If a source region is small, then the moments are likely noisy, so we zero
  // them. This is reasonable, given that small regions can get by with a flat
  // source approximation anyhow.
  if (source_regions_.is_small(sr)) {
    source_regions_.flux_moments_new(sr, g) = {0.0, 0.0, 0.0};
  } else {
    source_regions_.flux_moments_new(sr, g) *= (1.0 / volume);
  }
}

void LinearSourceDomain::set_flux_to_old_flux(int64_t sr, int g)
{
  source_regions_.scalar_flux_new(sr, g) =
    source_regions_.scalar_flux_old(sr, g);
  source_regions_.flux_moments_new(sr, g) =
    source_regions_.flux_moments_old(sr, g);
}

void LinearSourceDomain::accumulate_iteration_flux()
{
  // Accumulate scalar flux
  FlatSourceDomain::accumulate_iteration_flux();

  // Accumulate scalar flux moments
#pragma omp parallel for
  for (int64_t se = 0; se < n_source_elements(); se++) {
    source_regions_.flux_moments_t(se) += source_regions_.flux_moments_new(se);
  }
}

double LinearSourceDomain::evaluate_flux_at_point(
  Position r, int64_t sr, int g) const
{
  double phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g);

  Position local_r = r - source_regions_.centroid(sr);
  MomentArray phi_linear = source_regions_.flux_moments_t(sr, g);
  phi_linear *= 1.0 / (settings::n_batches - settings::n_inactive);

  MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();
  MomentArray phi_solved = invM * phi_linear;

  return phi_flat + phi_solved.dot(local_r);
}

} // namespace openmc

1	#include "openmc/random_ray/linear_source_domain.h"
2
3	#include "openmc/cell.h"
4	#include "openmc/geometry.h"
5	#include "openmc/material.h"
6	#include "openmc/message_passing.h"
7	#include "openmc/mgxs_interface.h"
8	#include "openmc/output.h"
9	#include "openmc/plot.h"
10	#include "openmc/random_ray/random_ray.h"
11	#include "openmc/simulation.h"
12	#include "openmc/tallies/filter.h"
13	#include "openmc/tallies/tally.h"
14	#include "openmc/tallies/tally_scoring.h"
15	#include "openmc/timer.h"
16
17	namespace openmc {
18
19	//==============================================================================
20	// LinearSourceDomain implementation
21	//==============================================================================
22
23	void LinearSourceDomain::batch_reset()	6,850✔
24	{
25	FlatSourceDomain::batch_reset();	6,850✔
26	#pragma omp parallel for	4,110✔
27	for (int64_t sr = 0; sr < n_source_regions(); sr++) {	10,483,272✔
28	source_regions_.centroid_iteration(sr) = {0.0, 0.0, 0.0};	10,480,532✔
29	source_regions_.mom_matrix(sr) = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};	10,480,532✔
30	}
31	#pragma omp parallel for	4,110✔
32	for (int64_t se = 0; se < n_source_elements(); se++) {	12,054,600✔
33	source_regions_.flux_moments_new(se) = {0.0, 0.0, 0.0};	12,051,860✔
34	}
35	}	6,850✔
36
37	void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh)	26,948,470✔
38	{
39	// Reset all source regions to zero (important for void regions)
40	for (int g = 0; g < negroups_; g++) {	57,888,140✔
41	srh.source(g) = 0.0;	30,939,670✔
42	}
43
44	// Add scattering + fission source
45	int material = srh.material();	26,948,470✔
46	int temp = srh.temperature_idx();	26,948,470✔
47	double density_mult = srh.density_mult();	26,948,470✔
48	if (material != MATERIAL_VOID) {	26,948,470✔
49	double inverse_k_eff = 1.0 / k_eff_;	26,548,470✔
50	MomentMatrix invM = srh.mom_matrix().inverse();	26,548,470✔
51
52	for (int g_out = 0; g_out < negroups_; g_out++) {	57,088,140✔
53	double sigma_t =
54	sigma_t_[(material * ntemperature_ + temp) * negroups_ + g_out] *	30,539,670✔
55	density_mult;	30,539,670✔
56
57	double scatter_flat = 0.0f;	30,539,670✔
58	double fission_flat = 0.0f;	30,539,670✔
59	MomentArray scatter_linear = {0.0, 0.0, 0.0};	30,539,670✔
60	MomentArray fission_linear = {0.0, 0.0, 0.0};	30,539,670✔
61
62	for (int g_in = 0; g_in < negroups_; g_in++) {	99,430,140✔
63	// Handles for the flat and linear components of the flux
64	double flux_flat = srh.scalar_flux_old(g_in);	68,890,470✔
65	MomentArray flux_linear = srh.flux_moments_old(g_in);	68,890,470✔
66
67	// Handles for cross sections
68	double sigma_s =
69	sigma_s_[((material * ntemperature_ + temp) * negroups_ + g_out) *	68,890,470✔
70	negroups_ +	68,890,470✔
71	g_in] *	68,890,470✔
72	density_mult;	68,890,470✔
73	double nu_sigma_f =
74	nu_sigma_f_[(material * ntemperature_ + temp) * negroups_ + g_in] *	68,890,470✔
75	density_mult;	68,890,470✔
76	double chi =
77	chi_[(material * ntemperature_ + temp) * negroups_ + g_out];	68,890,470✔
78
79	// Compute source terms for flat and linear components of the flux
80	scatter_flat += sigma_s * flux_flat;	68,890,470✔
81	scatter_linear += sigma_s * flux_linear;	68,890,470✔
82	if (settings::create_fission_neutrons) {	68,890,470!
83	fission_flat += nu_sigma_f * flux_flat * chi;	68,890,470✔
84	fission_linear += nu_sigma_f * flux_linear * chi;	68,890,470✔
85	}
86	}
87
88	// Compute the flat source term
89	srh.source(g_out) =	61,079,340✔
90	(scatter_flat + fission_flat * inverse_k_eff) / sigma_t;	30,539,670✔
91
92	// Compute the linear source terms. In the first 10 iterations when the
93	// centroids and spatial moments are not well known, we will leave the
94	// source gradients as zero so as to avoid causing any numerical
95	// instability. If a negative source is encountered, this region must be
96	// very small/noisy or have poorly developed spatial moments, so we zero
97	// the source gradients (effectively making this a flat source region
98	// temporarily), so as to improve stability.
99	if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) {	30,539,670✔
100	srh.source_gradients(g_out) =	18,690,020✔
101	invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);	37,380,040✔
102	} else {
103	srh.source_gradients(g_out) = {0.0, 0.0, 0.0};	11,849,650✔
104	}
105	}
106	}
107
108	// Add external source if in fixed source mode
109	if (settings::run_mode == RunMode::FIXED_SOURCE) {	26,948,470✔
110	for (int g = 0; g < negroups_; g++) {	56,080,540✔
111	srh.source(g) += srh.external_source(g);	29,345,270✔
112	}
113	}
114	}	26,948,470✔
115
116	void LinearSourceDomain::normalize_scalar_flux_and_volumes(	6,850✔
117	double total_active_distance_per_iteration)
118	{
119	double normalization_factor = 1.0 / total_active_distance_per_iteration;	6,850✔
120	double volume_normalization_factor =	6,850✔
121	1.0 / (total_active_distance_per_iteration * simulation::current_batch);	6,850✔
122
123	// Normalize flux to total distance travelled by all rays this iteration
124	#pragma omp parallel for	4,110✔
125	for (int64_t se = 0; se < n_source_elements(); se++) {	12,366,552✔
126	source_regions_.scalar_flux_new(se) *= normalization_factor;	12,363,812✔
127	source_regions_.flux_moments_new(se) *= normalization_factor;	12,363,812✔
128	}
129
130	// Accumulate cell-wise ray length tallies collected this iteration, then
131	// update the simulation-averaged cell-wise volume estimates
132	#pragma omp parallel for	4,110✔
133	for (int64_t sr = 0; sr < n_source_regions(); sr++) {	10,770,072✔
134	source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr);	10,767,332✔
135	source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr);	10,767,332✔
136	source_regions_.volume_t(sr) += source_regions_.volume(sr);	10,767,332✔
137	source_regions_.volume_sq_t(sr) += source_regions_.volume_sq(sr);	10,767,332✔
138	source_regions_.volume_naive(sr) =	21,534,664✔
139	source_regions_.volume(sr) * normalization_factor;	10,767,332✔
140	source_regions_.volume(sr) =	21,534,664✔
141	source_regions_.volume_t(sr) * volume_normalization_factor;	10,767,332✔
142	source_regions_.volume_sq(sr) =	21,534,664✔
143	source_regions_.volume_sq_t(sr) / source_regions_.volume_t(sr);	10,767,332✔
144	if (source_regions_.volume_t(sr) > 0.0) {	10,767,332✔
145	double inv_volume = 1.0 / source_regions_.volume_t(sr);	10,767,328✔
146	source_regions_.centroid(sr) = source_regions_.centroid_t(sr);	10,767,328✔
147	source_regions_.centroid(sr) *= inv_volume;	10,767,328✔
148	source_regions_.mom_matrix(sr) = source_regions_.mom_matrix_t(sr);	10,767,328✔
149	source_regions_.mom_matrix(sr) *= inv_volume;	10,767,328✔
150	}
151	}
152	}	6,850✔
153
154	void LinearSourceDomain::set_flux_to_flux_plus_source(	28,277,990✔
155	int64_t sr, double volume, int g)
156	{
157	int material = source_regions_.material(sr);	28,277,990✔
158	if (material == MATERIAL_VOID) {	28,277,990✔
159	FlatSourceDomain::set_flux_to_flux_plus_source(sr, volume, g);	400,000✔
160	} else {
161	source_regions_.scalar_flux_new(sr, g) /= volume;	27,877,990✔
162	source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);	27,877,990✔
163	}
164	// If a source region is small, then the moments are likely noisy, so we zero
165	// them. This is reasonable, given that small regions can get by with a flat
166	// source approximation anyhow.
167	if (source_regions_.is_small(sr)) {	28,277,990✔
168	source_regions_.flux_moments_new(sr, g) = {0.0, 0.0, 0.0};	2,871,140✔
169	} else {
170	source_regions_.flux_moments_new(sr, g) *= (1.0 / volume);	25,406,850✔
171	}
172	}	28,277,990✔
173
174	void LinearSourceDomain::set_flux_to_old_flux(int64_t sr, int g)	1,240✔
175	{
176	source_regions_.scalar_flux_new(sr, g) =	2,480✔
177	source_regions_.scalar_flux_old(sr, g);	1,240✔
178	source_regions_.flux_moments_new(sr, g) =	2,480✔
179	source_regions_.flux_moments_old(sr, g);	1,240✔
180	}	1,240✔
181
182	void LinearSourceDomain::accumulate_iteration_flux()	2,900✔
183	{
184	// Accumulate scalar flux
185	FlatSourceDomain::accumulate_iteration_flux();	2,900✔
186
187	// Accumulate scalar flux moments
188	#pragma omp parallel for	1,740✔
189	for (int64_t se = 0; se < n_source_elements(); se++) {	5,272,440✔
190	source_regions_.flux_moments_t(se) += source_regions_.flux_moments_new(se);	5,271,280✔
191	}
192	}	2,900✔
193
194	double LinearSourceDomain::evaluate_flux_at_point(	×
195	Position r, int64_t sr, int g) const
196	{
197	double phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g);	×
198
199	Position local_r = r - source_regions_.centroid(sr);	×
UNCOV 200	MomentArray phi_linear = source_regions_.flux_moments_t(sr, g);	×
UNCOV 201	phi_linear *= 1.0 / (settings::n_batches - settings::n_inactive);	×
202
UNCOV 203	MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();	×
UNCOV 204	MomentArray phi_solved = invM * phi_linear;	×
205
UNCOV 206	return phi_flat + phi_solved.dot(local_r);	×
207	}
208
209	} // namespace openmc

openmc-dev / openmc / 21887062786

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