15371300071

Committed 01 Jun 2025 05:04AM UTC coverage: 85.143% (+0.3%) from 84.827%

Build # 15371300071

Build Type

Pull #3176

github

Committed by

web-flow

Commit Message

Merge 4f739184a into cb95c784b

Pull Request Pull Request #3176: MeshFilter rotation - solution to issue #3166

Run Details

86 of 99 new or added lines in 4 files covered. (86.87%)

3707 existing lines in 117 files now uncovered.

52212 of 61323 relevant lines covered (85.14%)

42831974.38 hits per line

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

92.59

/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_neutron_source(double k_eff)
{
  simulation::time_update_src.start();

  double inverse_k_eff = 1.0 / k_eff;

// Reset all source regions to zero (important for void regions)
#pragma omp parallel for
  for (int64_t se = 0; se < n_source_elements(); se++) {
    source_regions_.source(se) = 0.0;
  }

#pragma omp parallel for
  for (int64_t sr = 0; sr < n_source_regions(); sr++) {
    int material = source_regions_.material(sr);
    if (material == MATERIAL_VOID) {
      continue;
    }
    MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();

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

      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 = source_regions_.scalar_flux_old(sr, g_in);
        MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in);

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

        // Compute source terms for flat and linear components of the flux
        scatter_flat += sigma_s * flux_flat;
        fission_flat += nu_sigma_f * flux_flat * chi;
        scatter_linear += sigma_s * flux_linear;
        fission_linear += nu_sigma_f * flux_linear * chi;
      }

      // Compute the flat source term
      source_regions_.source(sr, 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 &&
          source_regions_.source(sr, g_out) >= 0.0) {
        source_regions_.source_gradients(sr, g_out) =
          invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);
      } else {
        source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0};
      }
    }
  }

  if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Add external source to flat source term if in fixed source mode
#pragma omp parallel for
    for (int64_t se = 0; se < n_source_elements(); se++) {
      source_regions_.source(se) += source_regions_.external_source(se);
    }
  }

  simulation::time_update_src.stop();
}

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()	8,775✔
24	{
25	FlatSourceDomain::batch_reset();	8,775✔
26	#pragma omp parallel for	3,510✔
27	for (int64_t sr = 0; sr < n_source_regions(); sr++) {	23,708,790✔
28	source_regions_.centroid_iteration(sr) = {0.0, 0.0, 0.0};	23,703,525✔
29	source_regions_.mom_matrix(sr) = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};	23,703,525✔
30	}
31	#pragma omp parallel for	3,510✔
32	for (int64_t se = 0; se < n_source_elements(); se++) {	27,258,318✔
33	source_regions_.flux_moments_new(se) = {0.0, 0.0, 0.0};	27,253,053✔
34	}
35	}	8,775✔
36
37	void LinearSourceDomain::update_neutron_source(double k_eff)	8,775✔
38	{
39	simulation::time_update_src.start();	8,775✔
40
41	double inverse_k_eff = 1.0 / k_eff;	8,775✔
42
43	// Reset all source regions to zero (important for void regions)
44	#pragma omp parallel for	3,510✔
45	for (int64_t se = 0; se < n_source_elements(); se++) {	27,258,318✔
46	source_regions_.source(se) = 0.0;	27,253,053✔
47	}
48
49	#pragma omp parallel for	3,510✔
50	for (int64_t sr = 0; sr < n_source_regions(); sr++) {	23,708,790✔
51	int material = source_regions_.material(sr);	23,703,525✔
52	if (material == MATERIAL_VOID) {	23,703,525✔
53	continue;	360,000✔
54	}
55	MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();	23,343,525✔
56
57	for (int g_out = 0; g_out < negroups_; g_out++) {	50,236,578✔
58	double sigma_t = sigma_t_[material * negroups_ + g_out];	26,893,053✔
59
60	double scatter_flat = 0.0f;	26,893,053✔
61	double fission_flat = 0.0f;	26,893,053✔
62	MomentArray scatter_linear = {0.0, 0.0, 0.0};	26,893,053✔
63	MomentArray fission_linear = {0.0, 0.0, 0.0};	26,893,053✔
64
65	for (int g_in = 0; g_in < negroups_; g_in++) {	87,655,230✔
66	// Handles for the flat and linear components of the flux
67	double flux_flat = source_regions_.scalar_flux_old(sr, g_in);	60,762,177✔
68	MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in);	60,762,177✔
69
70	// Handles for cross sections
71	double sigma_s =
72	sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in];	60,762,177✔
73	double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in];	60,762,177✔
74	double chi = chi_[material * negroups_ + g_out];	60,762,177✔
75
76	// Compute source terms for flat and linear components of the flux
77	scatter_flat += sigma_s * flux_flat;	60,762,177✔
78	fission_flat += nu_sigma_f * flux_flat * chi;	60,762,177✔
79	scatter_linear += sigma_s * flux_linear;	60,762,177✔
80	fission_linear += nu_sigma_f * flux_linear * chi;	60,762,177✔
81	}
82
83	// Compute the flat source term
84	source_regions_.source(sr, g_out) =	53,786,106✔
85	(scatter_flat + fission_flat * inverse_k_eff) / sigma_t;	26,893,053✔
86
87	// Compute the linear source terms. In the first 10 iterations when the
88	// centroids and spatial moments are not well known, we will leave the
89	// source gradients as zero so as to avoid causing any numerical
90	// instability. If a negative source is encountered, this region must be
91	// very small/noisy or have poorly developed spatial moments, so we zero
92	// the source gradients (effectively making this a flat source region
93	// temporarily), so as to improve stability.
94	if (simulation::current_batch > 10 &&	46,288,953✔
95	source_regions_.source(sr, g_out) >= 0.0) {	19,395,900✔
96	source_regions_.source_gradients(sr, g_out) =	16,821,018✔
97	invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t);	33,642,036✔
98	} else {
99	source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0};	10,072,035✔
100	}
101	}
102	}
103
104	if (settings::run_mode == RunMode::FIXED_SOURCE) {	8,775✔
105	// Add external source to flat source term if in fixed source mode
106	#pragma omp parallel for	2,730✔
107	for (int64_t se = 0; se < n_source_elements(); se++) {	25,875,864✔
108	source_regions_.source(se) += source_regions_.external_source(se);	25,871,769✔
109	}
110	}
111
112	simulation::time_update_src.stop();	8,775✔
113	}	8,775✔
114
115	void LinearSourceDomain::normalize_scalar_flux_and_volumes(	8,775✔
116	double total_active_distance_per_iteration)
117	{
118	double normalization_factor = 1.0 / total_active_distance_per_iteration;	8,775✔
119	double volume_normalization_factor =	8,775✔
120	1.0 / (total_active_distance_per_iteration * simulation::current_batch);	8,775✔
121
122	// Normalize flux to total distance travelled by all rays this iteration
123	#pragma omp parallel for	3,510✔
124	for (int64_t se = 0; se < n_source_elements(); se++) {	27,776,322✔
125	source_regions_.scalar_flux_new(se) *= normalization_factor;	27,771,057✔
126	source_regions_.flux_moments_new(se) *= normalization_factor;	27,771,057✔
127	}
128
129	// Accumulate cell-wise ray length tallies collected this iteration, then
130	// update the simulation-averaged cell-wise volume estimates
131	#pragma omp parallel for	3,510✔
132	for (int64_t sr = 0; sr < n_source_regions(); sr++) {	24,220,962✔
133	source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr);	24,215,697✔
134	source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr);	24,215,697✔
135	source_regions_.volume_t(sr) += source_regions_.volume(sr);	24,215,697✔
136	source_regions_.volume_sq_t(sr) += source_regions_.volume_sq(sr);	24,215,697✔
137	source_regions_.volume_naive(sr) =	48,431,394✔
138	source_regions_.volume(sr) * normalization_factor;	24,215,697✔
139	source_regions_.volume(sr) =	48,431,394✔
140	source_regions_.volume_t(sr) * volume_normalization_factor;	24,215,697✔
141	source_regions_.volume_sq(sr) =	48,431,394✔
142	source_regions_.volume_sq_t(sr) / source_regions_.volume_t(sr);	24,215,697✔
143	if (source_regions_.volume_t(sr) > 0.0) {	24,215,697✔
144	double inv_volume = 1.0 / source_regions_.volume_t(sr);	24,215,688✔
145	source_regions_.centroid(sr) = source_regions_.centroid_t(sr);	24,215,688✔
146	source_regions_.centroid(sr) *= inv_volume;	24,215,688✔
147	source_regions_.mom_matrix(sr) = source_regions_.mom_matrix_t(sr);	24,215,688✔
148	source_regions_.mom_matrix(sr) *= inv_volume;	24,215,688✔
149	}
150	}
151	}	8,775✔
152
153	void LinearSourceDomain::set_flux_to_flux_plus_source(	42,337,785✔
154	int64_t sr, double volume, int g)
155	{
156	int material = source_regions_.material(sr);	42,337,785✔
157	if (material == MATERIAL_VOID) {	42,337,785✔
158	FlatSourceDomain::set_flux_to_flux_plus_source(sr, volume, g);	600,000✔
159	} else {
160	source_regions_.scalar_flux_new(sr, g) /= volume;	41,737,785✔
161	source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g);	41,737,785✔
162	}
163	// If a source region is small, then the moments are likely noisy, so we zero
164	// them. This is reasonable, given that small regions can get by with a flat
165	// source approximation anyhow.
166	if (source_regions_.is_small(sr)) {	42,337,785✔
167	source_regions_.flux_moments_new(sr, g) = {0.0, 0.0, 0.0};	4,306,710✔
168	} else {
169	source_regions_.flux_moments_new(sr, g) *= (1.0 / volume);	38,031,075✔
170	}
171	}	42,337,785✔
172
173	void LinearSourceDomain::set_flux_to_old_flux(int64_t sr, int g)	1,350✔
174	{
175	source_regions_.scalar_flux_new(sr, g) =	2,700✔
176	source_regions_.scalar_flux_old(sr, g);	1,350✔
177	source_regions_.flux_moments_new(sr, g) =	2,700✔
178	source_regions_.flux_moments_old(sr, g);	1,350✔
179	}	1,350✔
180
181	void LinearSourceDomain::accumulate_iteration_flux()	3,600✔
182	{
183	// Accumulate scalar flux
184	FlatSourceDomain::accumulate_iteration_flux();	3,600✔
185
186	// Accumulate scalar flux moments
187	#pragma omp parallel for	1,440✔
188	for (int64_t se = 0; se < n_source_elements(); se++) {	11,838,780✔
189	source_regions_.flux_moments_t(se) += source_regions_.flux_moments_new(se);	11,836,620✔
190	}
191	}	3,600✔
192
UNCOV 193	double LinearSourceDomain::evaluate_flux_at_point(	×
194	Position r, int64_t sr, int g) const
195	{
UNCOV 196	double phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g);	×
197
UNCOV 198	Position local_r = r - source_regions_.centroid(sr);	×
UNCOV 199	MomentArray phi_linear = source_regions_.flux_moments_t(sr, g);	×
UNCOV 200	phi_linear *= 1.0 / (settings::n_batches - settings::n_inactive);	×
201
UNCOV 202	MomentMatrix invM = source_regions_.mom_matrix(sr).inverse();	×
UNCOV 203	MomentArray phi_solved = invM * phi_linear;	×
204
UNCOV 205	return phi_flat + phi_solved.dot(local_r);	×
206	}
207
208	} // namespace openmc

openmc-dev / openmc / 15371300071

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