• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 21458205742

28 Jan 2026 10:40PM UTC coverage: 81.974% (-0.02%) from 81.994%
21458205742

Pull #3751

github

web-flow
Merge a14249281 into 008d58460
Pull Request #3751: Resolve conflict with weight windows and global russian roulette

17239 of 24006 branches covered (71.81%)

Branch coverage included in aggregate %.

87 of 100 new or added lines in 5 files covered. (87.0%)

94 existing lines in 2 files now uncovered.

55702 of 64975 relevant lines covered (85.73%)

44017310.75 hits per line

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

78.66
/src/weight_windows.cpp
1
#include "openmc/weight_windows.h"
2

3
#include <algorithm>
4
#include <cassert>
5
#include <cmath>
6
#include <set>
7
#include <string>
8

9
#include "xtensor/xdynamic_view.hpp"
10
#include "xtensor/xindex_view.hpp"
11
#include "xtensor/xio.hpp"
12
#include "xtensor/xmasked_view.hpp"
13
#include "xtensor/xnoalias.hpp"
14
#include "xtensor/xview.hpp"
15

16
#include "openmc/error.h"
17
#include "openmc/file_utils.h"
18
#include "openmc/hdf5_interface.h"
19
#include "openmc/mesh.h"
20
#include "openmc/message_passing.h"
21
#include "openmc/nuclide.h"
22
#include "openmc/output.h"
23
#include "openmc/particle.h"
24
#include "openmc/particle_data.h"
25
#include "openmc/physics_common.h"
26
#include "openmc/random_ray/flat_source_domain.h"
27
#include "openmc/search.h"
28
#include "openmc/settings.h"
29
#include "openmc/simulation.h"
30
#include "openmc/tallies/filter_energy.h"
31
#include "openmc/tallies/filter_mesh.h"
32
#include "openmc/tallies/filter_particle.h"
33
#include "openmc/tallies/tally.h"
34
#include "openmc/xml_interface.h"
35

36
#include <fmt/core.h>
37

38
namespace openmc {
39

40
//==============================================================================
41
// Global variables
42
//==============================================================================
43

44
namespace variance_reduction {
45

46
std::unordered_map<int32_t, int32_t> ww_map;
47
openmc::vector<unique_ptr<WeightWindows>> weight_windows;
48
openmc::vector<unique_ptr<WeightWindowsGenerator>> weight_windows_generators;
49

50
} // namespace variance_reduction
51

52
//==============================================================================
53
// WeightWindowSettings implementation
54
//==============================================================================
55

56
WeightWindows::WeightWindows(int32_t id)
247✔
57
{
58
  index_ = variance_reduction::weight_windows.size();
247✔
59
  set_id(id);
247✔
60
  set_defaults();
247✔
61
}
247✔
62

63
WeightWindows::WeightWindows(pugi::xml_node node)
92✔
64
{
65
  // Make sure required elements are present
66
  const vector<std::string> required_elems {
67
    "id", "particle_type", "lower_ww_bounds", "upper_ww_bounds"};
644✔
68
  for (const auto& elem : required_elems) {
460✔
69
    if (!check_for_node(node, elem.c_str())) {
368!
70
      fatal_error(fmt::format("Must specify <{}> for weight windows.", elem));
×
71
    }
72
  }
73

74
  // Get weight windows ID
75
  int32_t id = std::stoi(get_node_value(node, "id"));
92✔
76
  this->set_id(id);
92✔
77

78
  // get the particle type
79
  auto particle_type_str = std::string(get_node_value(node, "particle_type"));
92✔
80
  particle_type_ = openmc::str_to_particle_type(particle_type_str);
92✔
81

82
  // Determine associated mesh
83
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
92✔
84
  set_mesh(model::mesh_map.at(mesh_id));
92✔
85

86
  // energy bounds
87
  if (check_for_node(node, "energy_bounds"))
92✔
88
    energy_bounds_ = get_node_array<double>(node, "energy_bounds");
77✔
89

90
  // get the survival value - optional
91
  if (check_for_node(node, "survival_ratio")) {
92!
92
    survival_ratio_ = std::stod(get_node_value(node, "survival_ratio"));
92✔
93
    if (survival_ratio_ <= 1)
92!
94
      fatal_error("Survival to lower weight window ratio must bigger than 1 "
×
95
                  "and less than the upper to lower weight window ratio.");
96
  }
97

98
  // get the max lower bound ratio - optional
99
  if (check_for_node(node, "max_lower_bound_ratio")) {
92✔
100
    max_lb_ratio_ = std::stod(get_node_value(node, "max_lower_bound_ratio"));
34✔
101
    if (max_lb_ratio_ < 1.0) {
34!
102
      fatal_error("Maximum lower bound ratio must be larger than 1");
×
103
    }
104
  }
105

106
  // get the max split - optional
107
  if (check_for_node(node, "max_split")) {
92!
108
    max_split_ = std::stod(get_node_value(node, "max_split"));
92✔
109
    if (max_split_ <= 1)
92!
110
      fatal_error("max split must be larger than 1");
×
111
  }
112

113
  // weight cutoff - optional
114
  if (check_for_node(node, "weight_cutoff")) {
92!
115
    weight_cutoff_ = std::stod(get_node_value(node, "weight_cutoff"));
92✔
116
    if (weight_cutoff_ <= 0)
92!
117
      fatal_error("weight_cutoff must be larger than 0");
×
118
    if (weight_cutoff_ > 1)
92!
119
      fatal_error("weight_cutoff must be less than 1");
×
120
  }
121

122
  // read the lower/upper weight bounds
123
  this->set_bounds(get_node_array<double>(node, "lower_ww_bounds"),
92✔
124
    get_node_array<double>(node, "upper_ww_bounds"));
184✔
125

126
  set_defaults();
92✔
127
}
92✔
128

129
WeightWindows::~WeightWindows()
339✔
130
{
131
  variance_reduction::ww_map.erase(id());
339✔
132
}
339✔
133

134
WeightWindows* WeightWindows::create(int32_t id)
93✔
135
{
136
  variance_reduction::weight_windows.push_back(make_unique<WeightWindows>());
93✔
137
  auto wws = variance_reduction::weight_windows.back().get();
93✔
138
  variance_reduction::ww_map[wws->id()] =
93✔
139
    variance_reduction::weight_windows.size() - 1;
93✔
140
  return wws;
93✔
141
}
142

143
WeightWindows* WeightWindows::from_hdf5(
11✔
144
  hid_t wws_group, const std::string& group_name)
145
{
146
  // collect ID from the name of this group
147
  hid_t ww_group = open_group(wws_group, group_name);
11✔
148

149
  auto wws = WeightWindows::create();
11✔
150

151
  std::string particle_type;
11✔
152
  read_dataset(ww_group, "particle_type", particle_type);
11✔
153
  wws->particle_type_ = openmc::str_to_particle_type(particle_type);
11✔
154

155
  read_dataset<double>(ww_group, "energy_bounds", wws->energy_bounds_);
11✔
156

157
  int32_t mesh_id;
158
  read_dataset(ww_group, "mesh", mesh_id);
11✔
159

160
  if (model::mesh_map.count(mesh_id) == 0) {
11!
161
    fatal_error(
×
162
      fmt::format("Mesh {} used in weight windows does not exist.", mesh_id));
×
163
  }
164
  wws->set_mesh(model::mesh_map[mesh_id]);
11✔
165

166
  wws->lower_ww_ = xt::empty<double>(wws->bounds_size());
11✔
167
  wws->upper_ww_ = xt::empty<double>(wws->bounds_size());
11✔
168

169
  read_dataset<double>(ww_group, "lower_ww_bounds", wws->lower_ww_);
11✔
170
  read_dataset<double>(ww_group, "upper_ww_bounds", wws->upper_ww_);
11✔
171
  read_dataset(ww_group, "survival_ratio", wws->survival_ratio_);
11✔
172
  read_dataset(ww_group, "max_lower_bound_ratio", wws->max_lb_ratio_);
11✔
173
  read_dataset(ww_group, "max_split", wws->max_split_);
11✔
174
  read_dataset(ww_group, "weight_cutoff", wws->weight_cutoff_);
11✔
175

176
  close_group(ww_group);
11✔
177

178
  return wws;
11✔
179
}
11✔
180

181
void WeightWindows::set_defaults()
421✔
182
{
183
  // set energy bounds to the min/max energy supported by the data
184
  if (energy_bounds_.size() == 0) {
421✔
185
    int p_type = static_cast<int>(particle_type_);
262✔
186
    energy_bounds_.push_back(data::energy_min[p_type]);
262✔
187
    energy_bounds_.push_back(data::energy_max[p_type]);
262✔
188
  }
189
}
421✔
190

191
void WeightWindows::allocate_ww_bounds()
553✔
192
{
193
  auto shape = bounds_size();
553✔
194
  if (shape[0] * shape[1] == 0) {
553!
195
    auto msg = fmt::format(
196
      "Size of weight window bounds is zero for WeightWindows {}", id());
×
197
    warning(msg);
×
198
  }
×
199
  lower_ww_ = xt::empty<double>(shape);
553✔
200
  lower_ww_.fill(-1);
553✔
201
  upper_ww_ = xt::empty<double>(shape);
553✔
202
  upper_ww_.fill(-1);
553✔
203
}
553✔
204

205
void WeightWindows::set_id(int32_t id)
493✔
206
{
207
  assert(id >= 0 || id == C_NONE);
402!
208

209
  // Clear entry in mesh map in case one was already assigned
210
  if (id_ != C_NONE) {
493!
211
    variance_reduction::ww_map.erase(id_);
493✔
212
    id_ = C_NONE;
493✔
213
  }
214

215
  // Ensure no other mesh has the same ID
216
  if (variance_reduction::ww_map.find(id) != variance_reduction::ww_map.end()) {
493!
217
    throw std::runtime_error {
×
218
      fmt::format("Two weight windows have the same ID: {}", id)};
×
219
  }
220

221
  // If no ID is specified, auto-assign the next ID in the sequence
222
  if (id == C_NONE) {
493✔
223
    id = 0;
247✔
224
    for (const auto& m : variance_reduction::weight_windows) {
269✔
225
      id = std::max(id, m->id_);
22✔
226
    }
227
    ++id;
247✔
228
  }
229

230
  // Update ID and entry in the mesh map
231
  id_ = id;
493✔
232
  variance_reduction::ww_map[id] = index_;
493✔
233
}
493✔
234

235
void WeightWindows::set_energy_bounds(span<const double> bounds)
214✔
236
{
237
  energy_bounds_.clear();
214✔
238
  energy_bounds_.insert(energy_bounds_.begin(), bounds.begin(), bounds.end());
214✔
239
  // if the mesh is set, allocate space for weight window bounds
240
  if (mesh_idx_ != C_NONE)
214!
241
    allocate_ww_bounds();
214✔
242
}
214✔
243

244
void WeightWindows::set_particle_type(ParticleType p_type)
258✔
245
{
246
  if (p_type != ParticleType::neutron && p_type != ParticleType::photon)
258!
247
    fatal_error(
×
248
      fmt::format("Particle type '{}' cannot be applied to weight windows.",
×
249
        particle_type_to_str(p_type)));
×
250
  particle_type_ = p_type;
258✔
251
}
258✔
252

253
void WeightWindows::set_mesh(int32_t mesh_idx)
339✔
254
{
255
  if (mesh_idx < 0 || mesh_idx >= model::meshes.size())
339!
256
    fatal_error(fmt::format("Could not find a mesh for index {}", mesh_idx));
×
257

258
  mesh_idx_ = mesh_idx;
339✔
259
  model::meshes[mesh_idx_]->prepare_for_point_location();
339✔
260
  allocate_ww_bounds();
339✔
261
}
339✔
262

263
void WeightWindows::set_mesh(const std::unique_ptr<Mesh>& mesh)
×
264
{
265
  set_mesh(mesh.get());
×
266
}
×
267

268
void WeightWindows::set_mesh(const Mesh* mesh)
×
269
{
270
  set_mesh(model::mesh_map[mesh->id_]);
×
271
}
×
272

273
const int WeightWindows::get_mesh_bin(const Particle& p) const
90,902,269✔
274
{
275
  // check for particle type
276
  if (particle_type_ != p.type())
90,902,269✔
277
    return C_NONE;
19,709,373✔
278

279
  // particle energy
280
  double E = p.E();
71,192,896✔
281

282
  // check to make sure energy is in range, expects sorted energy values
283
  if (E < energy_bounds_.front() || E > energy_bounds_.back())
71,192,896!
284
    return C_NONE;
92,818✔
285

286
  // Get mesh index for particle's position
287
  const auto& mesh = this->mesh();
71,100,078✔
288
  return mesh->get_bin(p.r());
71,100,078✔
289
}
290

NEW
291
WeightWindow WeightWindows::get_weight_window(const Particle& p) const
×
292
{
NEW
293
  return get_weight_window(p.E(), get_mesh_bin(p));
×
294
}
295

296
WeightWindow WeightWindows::get_weight_window(
71,085,244✔
297
  double E, const int mesh_bin) const
298
{
299
  // particle is outside the weight window mesh
300
  if (mesh_bin < 0)
71,085,244!
UNCOV
301
    return {};
×
302

303
  // get the mesh bin in energy group
304
  int energy_bin =
305
    lower_bound_index(energy_bounds_.begin(), energy_bounds_.end(), E);
71,085,244✔
306

307
  // mesh_bin += energy_bin * mesh->n_bins();
308
  // Create individual weight window
309
  WeightWindow ww;
71,085,244✔
310
  ww.lower_weight = lower_ww_(energy_bin, mesh_bin);
71,085,244✔
311
  ww.upper_weight = upper_ww_(energy_bin, mesh_bin);
71,085,244✔
312
  ww.survival_weight = ww.lower_weight * survival_ratio_;
71,085,244✔
313
  ww.max_lb_ratio = max_lb_ratio_;
71,085,244✔
314
  ww.max_split = max_split_;
71,085,244✔
315
  ww.weight_cutoff = weight_cutoff_;
71,085,244✔
316
  return ww;
71,085,244✔
317
}
318

319
std::array<int, 2> WeightWindows::bounds_size() const
781✔
320
{
321
  int num_spatial_bins = this->mesh()->n_bins();
781✔
322
  int num_energy_bins =
323
    energy_bounds_.size() > 0 ? energy_bounds_.size() - 1 : 1;
781✔
324
  return {num_energy_bins, num_spatial_bins};
781✔
325
}
326

327
template<class T>
328
void WeightWindows::check_bounds(const T& lower, const T& upper) const
103✔
329
{
330
  // make sure that the upper and lower bounds have the same size
331
  if (lower.size() != upper.size()) {
103!
332
    auto msg = fmt::format("The upper and lower weight window lengths do not "
×
333
                           "match.\n Lower size: {}\n Upper size: {}",
334
      lower.size(), upper.size());
×
335
    fatal_error(msg);
×
336
  }
×
337
  this->check_bounds(lower);
103✔
338
}
103✔
339

340
template<class T>
341
void WeightWindows::check_bounds(const T& bounds) const
103✔
342
{
343
  // check that the number of weight window entries is correct
344
  auto dims = this->bounds_size();
103✔
345
  if (bounds.size() != dims[0] * dims[1]) {
103!
346
    auto err_msg =
×
347
      fmt::format("In weight window domain {} the number of spatial "
348
                  "energy/spatial bins ({}) does not match the number "
349
                  "of weight bins ({})",
350
        id_, dims, bounds.size());
×
351
    fatal_error(err_msg);
×
352
  }
×
353
}
103✔
354

355
void WeightWindows::set_bounds(const xt::xtensor<double, 2>& lower_bounds,
×
356
  const xt::xtensor<double, 2>& upper_bounds)
357
{
358

359
  this->check_bounds(lower_bounds, upper_bounds);
×
360

361
  // set new weight window values
362
  lower_ww_ = lower_bounds;
×
363
  upper_ww_ = upper_bounds;
×
364
}
×
365

366
void WeightWindows::set_bounds(
×
367
  const xt::xtensor<double, 2>& lower_bounds, double ratio)
368
{
369
  this->check_bounds(lower_bounds);
×
370

371
  // set new weight window values
372
  lower_ww_ = lower_bounds;
×
373
  upper_ww_ = lower_bounds;
×
374
  upper_ww_ *= ratio;
×
375
}
×
376

377
void WeightWindows::set_bounds(
103✔
378
  span<const double> lower_bounds, span<const double> upper_bounds)
379
{
380
  check_bounds(lower_bounds, upper_bounds);
103✔
381
  auto shape = this->bounds_size();
103✔
382
  lower_ww_ = xt::empty<double>(shape);
103✔
383
  upper_ww_ = xt::empty<double>(shape);
103✔
384

385
  // set new weight window values
386
  xt::view(lower_ww_, xt::all()) =
206✔
387
    xt::adapt(lower_bounds.data(), lower_ww_.shape());
309✔
388
  xt::view(upper_ww_, xt::all()) =
206✔
389
    xt::adapt(upper_bounds.data(), upper_ww_.shape());
309✔
390
}
103✔
391

392
void WeightWindows::set_bounds(span<const double> lower_bounds, double ratio)
×
393
{
394
  this->check_bounds(lower_bounds);
×
395

396
  auto shape = this->bounds_size();
×
397
  lower_ww_ = xt::empty<double>(shape);
×
398
  upper_ww_ = xt::empty<double>(shape);
×
399

400
  // set new weight window values
401
  xt::view(lower_ww_, xt::all()) =
×
402
    xt::adapt(lower_bounds.data(), lower_ww_.shape());
×
403
  xt::view(upper_ww_, xt::all()) =
×
404
    xt::adapt(lower_bounds.data(), upper_ww_.shape());
×
405
  upper_ww_ *= ratio;
×
406
}
×
407

408
void WeightWindows::update_weights(const Tally* tally, const std::string& value,
252✔
409
  double threshold, double ratio, WeightWindowUpdateMethod method)
410
{
411
  ///////////////////////////
412
  // Setup and checks
413
  ///////////////////////////
414
  this->check_tally_update_compatibility(tally);
252✔
415

416
  // Dimensions of weight window arrays
417
  int e_bins = lower_ww_.shape()[0];
252✔
418
  int64_t mesh_bins = lower_ww_.shape()[1];
252✔
419

420
  // Initialize weight window arrays to -1.0 by default
421
#pragma omp parallel for collapse(2) schedule(static)
140✔
422
  for (int e = 0; e < e_bins; e++) {
934✔
423
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
424
      lower_ww_(e, m) = -1.0;
1,189,249✔
425
      upper_ww_(e, m) = -1.0;
1,189,249✔
426
    }
427
  }
428

429
  // determine which value to use
430
  const std::set<std::string> allowed_values = {"mean", "rel_err"};
1,260✔
431
  if (allowed_values.count(value) == 0) {
252!
432
    fatal_error(fmt::format("Invalid value '{}' specified for weight window "
×
433
                            "generation. Must be one of: 'mean' or 'rel_err'",
434
      value));
435
  }
436

437
  // determine the index of the specified score
438
  int score_index = tally->score_index("flux");
252✔
439
  if (score_index == C_NONE) {
252!
440
    fatal_error(
×
441
      fmt::format("A 'flux' score required for weight window generation "
×
442
                  "is not present on tally {}.",
443
        tally->id()));
×
444
  }
445

446
  ///////////////////////////
447
  // Extract tally data
448
  //
449
  // At the end of this section, the mean and rel_err array
450
  // is a 2D view of tally data (n_e_groups, n_mesh_bins)
451
  //
452
  ///////////////////////////
453

454
  // build a shape for a view of the tally results, this will always be
455
  // dimension 5 (3 filter dimensions, 1 score dimension, 1 results dimension)
456
  // Look for the size of the last dimension of the results array
457
  const auto& results_arr = tally->results();
252✔
458
  const int results_dim = static_cast<int>(results_arr.shape()[2]);
252✔
459
  std::array<int, 5> shape = {1, 1, 1, tally->n_scores(), results_dim};
252✔
460

461
  // set the shape for the filters applied on the tally
462
  for (int i = 0; i < tally->filters().size(); i++) {
964✔
463
    const auto& filter = model::tally_filters[tally->filters(i)];
712✔
464
    shape[i] = filter->n_bins();
712✔
465
  }
466

467
  // build the transpose information to re-order data according to filter type
468
  std::array<int, 5> transpose = {0, 1, 2, 3, 4};
252✔
469

470
  // track our filter types and where we've added new ones
471
  std::vector<FilterType> filter_types = tally->filter_types();
252✔
472

473
  // assign other filter types to dummy positions if needed
474
  if (!tally->has_filter(FilterType::PARTICLE))
252✔
475
    filter_types.push_back(FilterType::PARTICLE);
22✔
476

477
  if (!tally->has_filter(FilterType::ENERGY))
252✔
478
    filter_types.push_back(FilterType::ENERGY);
22✔
479

480
  // particle axis mapping
481
  transpose[0] =
252✔
482
    std::find(filter_types.begin(), filter_types.end(), FilterType::PARTICLE) -
252✔
483
    filter_types.begin();
252✔
484

485
  // energy axis mapping
486
  transpose[1] =
252✔
487
    std::find(filter_types.begin(), filter_types.end(), FilterType::ENERGY) -
252✔
488
    filter_types.begin();
252✔
489

490
  // mesh axis mapping
491
  transpose[2] =
252✔
492
    std::find(filter_types.begin(), filter_types.end(), FilterType::MESH) -
252✔
493
    filter_types.begin();
252✔
494

495
  // get a fully reshaped view of the tally according to tally ordering of
496
  // filters
497
  auto tally_values = xt::reshape_view(results_arr, shape);
252✔
498

499
  // get a that is (particle, energy, mesh, scores, values)
500
  auto transposed_view = xt::transpose(tally_values, transpose);
252✔
501

502
  // determine the dimension and index of the particle data
503
  int particle_idx = 0;
252✔
504
  if (tally->has_filter(FilterType::PARTICLE)) {
252✔
505
    // get the particle filter
506
    auto pf = tally->get_filter<ParticleFilter>();
230✔
507
    const auto& particles = pf->particles();
230✔
508

509
    // find the index of the particle that matches these weight windows
510
    auto p_it =
511
      std::find(particles.begin(), particles.end(), this->particle_type_);
230✔
512
    // if the particle filter doesn't have particle data for the particle
513
    // used on this weight windows instance, report an error
514
    if (p_it == particles.end()) {
230!
515
      auto msg = fmt::format("Particle type '{}' not present on Filter {} for "
516
                             "Tally {} used to update WeightWindows {}",
517
        particle_type_to_str(this->particle_type_), pf->id(), tally->id(),
×
518
        this->id());
×
519
      fatal_error(msg);
×
520
    }
×
521

522
    // use the index of the particle in the filter to down-select data later
523
    particle_idx = p_it - particles.begin();
230✔
524
  }
525

526
  // down-select data based on particle and score
527
  auto sum = xt::dynamic_view(
1,260✔
528
    transposed_view, {particle_idx, xt::all(), xt::all(), score_index,
504✔
529
                       static_cast<int>(TallyResult::SUM)});
1,008✔
530
  auto sum_sq = xt::dynamic_view(
1,260✔
531
    transposed_view, {particle_idx, xt::all(), xt::all(), score_index,
504✔
532
                       static_cast<int>(TallyResult::SUM_SQ)});
1,008✔
533
  int n = tally->n_realizations_;
252✔
534

535
  //////////////////////////////////////////////
536
  //
537
  // Assign new weight windows
538
  //
539
  // Use references to the existing weight window data
540
  // to store and update the values
541
  //
542
  //////////////////////////////////////////////
543

544
  // up to this point the data arrays are views into the tally results (no
545
  // computation has been performed) now we'll switch references to the tally's
546
  // bounds to avoid allocating additional memory
547
  auto& new_bounds = this->lower_ww_;
252✔
548
  auto& rel_err = this->upper_ww_;
252✔
549

550
  // get mesh volumes
551
  auto mesh_vols = this->mesh()->volumes();
252✔
552

553
  // Calculate mean (new_bounds) and relative error
554
#pragma omp parallel for collapse(2) schedule(static)
140✔
555
  for (int e = 0; e < e_bins; e++) {
934✔
556
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
557
      // Calculate mean
558
      new_bounds(e, m) = sum(e, m) / n;
1,189,249✔
559
      // Calculate relative error
560
      if (sum(e, m) > 0.0) {
1,189,249✔
561
        double mean_val = new_bounds(e, m);
101,480✔
562
        double variance = (sum_sq(e, m) / n - mean_val * mean_val) / (n - 1);
101,480✔
563
        rel_err(e, m) = std::sqrt(variance) / mean_val;
101,480✔
564
      } else {
565
        rel_err(e, m) = INFTY;
1,087,769✔
566
      }
567
      if (value == "rel_err") {
1,189,249✔
568
        new_bounds(e, m) = 1.0 / rel_err(e, m);
345,000✔
569
      }
570
    }
571
  }
572

573
  // Divide by volume of mesh elements
574
#pragma omp parallel for collapse(2) schedule(static)
140✔
575
  for (int e = 0; e < e_bins; e++) {
934✔
576
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
577
      new_bounds(e, m) /= mesh_vols[m];
1,189,249✔
578
    }
579
  }
580

581
  if (method == WeightWindowUpdateMethod::MAGIC) {
252✔
582
    // For MAGIC, weight windows are proportional to the forward fluxes.
583
    // We normalize weight windows independently for each energy group.
584

585
    // Find group maximum and normalize (per energy group)
586
    for (int e = 0; e < e_bins; e++) {
1,870✔
587
      double group_max = 0.0;
1,716✔
588

589
      // Find maximum value across all elements in this energy group
590
#pragma omp parallel for schedule(static) reduction(max : group_max)
936✔
591
      for (int64_t m = 0; m < mesh_bins; m++) {
1,092,505✔
592
        if (new_bounds(e, m) > group_max) {
1,091,725✔
593
          group_max = new_bounds(e, m);
2,520✔
594
        }
595
      }
596

597
      // Normalize values in this energy group by the maximum value
598
      if (group_max > 0.0) {
1,716✔
599
        double norm_factor = 1.0 / (2.0 * group_max);
1,683✔
600
#pragma omp parallel for schedule(static)
918✔
601
        for (int64_t m = 0; m < mesh_bins; m++) {
1,091,590✔
602
          new_bounds(e, m) *= norm_factor;
1,090,825✔
603
        }
604
      }
605
    }
606
  } else {
607
    // For FW-CADIS, weight windows are inversely proportional to the adjoint
608
    // fluxes. We normalize the weight windows across all energy groups.
609
#pragma omp parallel for collapse(2) schedule(static)
56✔
610
    for (int e = 0; e < e_bins; e++) {
84✔
611
      for (int64_t m = 0; m < mesh_bins; m++) {
97,566✔
612
        // Take the inverse, but are careful not to divide by zero
613
        if (new_bounds(e, m) != 0.0) {
97,524✔
614
          new_bounds(e, m) = 1.0 / new_bounds(e, m);
69,660✔
615
        } else {
616
          new_bounds(e, m) = 0.0;
27,864!
617
        }
618
      }
619
    }
620

621
    // Find the maximum value across all elements
622
    double max_val = 0.0;
98✔
623
#pragma omp parallel for collapse(2) schedule(static) reduction(max : max_val)
56✔
624
    for (int e = 0; e < e_bins; e++) {
84✔
625
      for (int64_t m = 0; m < mesh_bins; m++) {
97,566✔
626
        if (new_bounds(e, m) > max_val) {
97,524✔
627
          max_val = new_bounds(e, m);
405✔
628
        }
629
      }
630
    }
631

632
    // Parallel normalization
633
    if (max_val > 0.0) {
98✔
634
      double norm_factor = 1.0 / (2.0 * max_val);
68✔
635
#pragma omp parallel for collapse(2) schedule(static)
38✔
636
      for (int e = 0; e < e_bins; e++) {
60✔
637
        for (int64_t m = 0; m < mesh_bins; m++) {
69,690✔
638
          new_bounds(e, m) *= norm_factor;
69,660✔
639
        }
640
      }
641
    }
642
  }
643

644
  // Final processing
645
#pragma omp parallel for collapse(2) schedule(static)
140✔
646
  for (int e = 0; e < e_bins; e++) {
934✔
647
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
648
      // Values where the mean is zero should be ignored
649
      if (sum(e, m) <= 0.0) {
1,189,249✔
650
        new_bounds(e, m) = -1.0;
1,087,769✔
651
      }
652
      // Values where the relative error is higher than the threshold should be
653
      // ignored
654
      else if (rel_err(e, m) > threshold) {
101,480✔
655
        new_bounds(e, m) = -1.0;
1,420✔
656
      }
657
      // Set the upper bounds
658
      upper_ww_(e, m) = ratio * lower_ww_(e, m);
1,189,249✔
659
    }
660
  }
661
}
252✔
662

663
void WeightWindows::check_tally_update_compatibility(const Tally* tally)
252✔
664
{
665
  // define the set of allowed filters for the tally
666
  const std::set<FilterType> allowed_filters = {
667
    FilterType::MESH, FilterType::ENERGY, FilterType::PARTICLE};
252✔
668

669
  // retrieve a mapping of filter type to filter index for the tally
670
  auto filter_indices = tally->filter_indices();
252✔
671

672
  // a mesh filter is required for a tally used to update weight windows
673
  if (!filter_indices.count(FilterType::MESH)) {
252!
674
    fatal_error(
×
675
      "A mesh filter is required for a tally to update weight window bounds");
676
  }
677

678
  // ensure the mesh filter is using the same mesh as this weight window object
679
  auto mesh_filter = tally->get_filter<MeshFilter>();
252✔
680

681
  // make sure that all of the filters present on the tally are allowed
682
  for (auto filter_pair : filter_indices) {
964✔
683
    if (allowed_filters.find(filter_pair.first) == allowed_filters.end()) {
712!
684
      fatal_error(fmt::format("Invalid filter type '{}' found on tally "
×
685
                              "used for weight window generation.",
686
        model::tally_filters[tally->filters(filter_pair.second)]->type_str()));
×
687
    }
688
  }
689

690
  if (mesh_filter->mesh() != mesh_idx_) {
252!
691
    int32_t mesh_filter_id = model::meshes[mesh_filter->mesh()]->id();
×
692
    int32_t ww_mesh_id = model::meshes[this->mesh_idx_]->id();
×
693
    fatal_error(fmt::format("Mesh filter {} uses a different mesh ({}) than "
×
694
                            "weight window {} mesh ({})",
695
      mesh_filter->id(), mesh_filter_id, id_, ww_mesh_id));
×
696
  }
697

698
  // if an energy filter exists, make sure the energy grid matches that of this
699
  // weight window object
700
  if (auto energy_filter = tally->get_filter<EnergyFilter>()) {
252✔
701
    std::vector<double> filter_bins = energy_filter->bins();
230✔
702
    std::set<double> filter_e_bounds(
703
      energy_filter->bins().begin(), energy_filter->bins().end());
230✔
704
    if (filter_e_bounds.size() != energy_bounds().size()) {
230!
705
      fatal_error(
×
706
        fmt::format("Energy filter {} does not have the same number of energy "
×
707
                    "bounds ({}) as weight window object {} ({})",
708
          energy_filter->id(), filter_e_bounds.size(), id_,
×
709
          energy_bounds().size()));
×
710
    }
711

712
    for (auto e : energy_bounds()) {
2,252✔
713
      if (filter_e_bounds.count(e) == 0) {
2,022!
714
        fatal_error(fmt::format(
×
715
          "Energy bounds of filter {} and weight windows {} do not match",
716
          energy_filter->id(), id_));
×
717
      }
718
    }
719
  }
230✔
720
}
252✔
721

722
void WeightWindows::to_hdf5(hid_t group) const
134✔
723
{
724
  hid_t ww_group = create_group(group, fmt::format("weight_windows_{}", id()));
268✔
725

726
  write_dataset(ww_group, "mesh", this->mesh()->id());
134✔
727
  write_dataset(
134✔
728
    ww_group, "particle_type", openmc::particle_type_to_str(particle_type_));
268✔
729
  write_dataset(ww_group, "energy_bounds", energy_bounds_);
134✔
730
  write_dataset(ww_group, "lower_ww_bounds", lower_ww_);
134✔
731
  write_dataset(ww_group, "upper_ww_bounds", upper_ww_);
134✔
732
  write_dataset(ww_group, "survival_ratio", survival_ratio_);
134✔
733
  write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_);
134✔
734
  write_dataset(ww_group, "max_split", max_split_);
134✔
735
  write_dataset(ww_group, "weight_cutoff", weight_cutoff_);
134✔
736

737
  close_group(ww_group);
134✔
738
}
134✔
739

740
WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
82✔
741
{
742
  // read information from the XML node
743
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
82✔
744
  int32_t mesh_idx = model::mesh_map[mesh_id];
82✔
745
  max_realizations_ = std::stoi(get_node_value(node, "max_realizations"));
82✔
746

747
  int32_t active_batches = settings::n_batches - settings::n_inactive;
82✔
748
  if (max_realizations_ > active_batches) {
82✔
749
    auto msg =
750
      fmt::format("The maximum number of specified tally realizations ({}) is "
751
                  "greater than the number of active batches ({}).",
752
        max_realizations_, active_batches);
31✔
753
    warning(msg);
17✔
754
  }
17✔
755
  auto tmp_str = get_node_value(node, "particle_type", true, true);
82✔
756
  auto particle_type = str_to_particle_type(tmp_str);
82✔
757

758
  update_interval_ = std::stoi(get_node_value(node, "update_interval"));
82✔
759
  on_the_fly_ = get_node_value_bool(node, "on_the_fly");
82✔
760

761
  std::vector<double> e_bounds;
82✔
762
  if (check_for_node(node, "energy_bounds")) {
82✔
763
    e_bounds = get_node_array<double>(node, "energy_bounds");
23✔
764
  } else {
765
    int p_type = static_cast<int>(particle_type);
59✔
766
    e_bounds.push_back(data::energy_min[p_type]);
59✔
767
    e_bounds.push_back(data::energy_max[p_type]);
59✔
768
  }
769

770
  // set method
771
  std::string method_string = get_node_value(node, "method");
82✔
772
  if (method_string == "magic") {
82✔
773
    method_ = WeightWindowUpdateMethod::MAGIC;
33✔
774
    if (settings::solver_type == SolverType::RANDOM_RAY &&
33!
775
        FlatSourceDomain::adjoint_) {
776
      fatal_error("Random ray weight window generation with MAGIC cannot be "
×
777
                  "done in adjoint mode.");
778
    }
779
  } else if (method_string == "fw_cadis") {
49!
780
    method_ = WeightWindowUpdateMethod::FW_CADIS;
49✔
781
    if (settings::solver_type != SolverType::RANDOM_RAY) {
49!
782
      fatal_error("FW-CADIS can only be run in random ray solver mode.");
×
783
    }
784
    FlatSourceDomain::adjoint_ = true;
49✔
785
  } else {
786
    fatal_error(fmt::format(
×
787
      "Unknown weight window update method '{}' specified", method_string));
788
  }
789

790
  // parse non-default update parameters if specified
791
  if (check_for_node(node, "update_parameters")) {
82✔
792
    pugi::xml_node params_node = node.child("update_parameters");
22✔
793
    if (check_for_node(params_node, "value"))
22!
794
      tally_value_ = get_node_value(params_node, "value");
22✔
795
    if (check_for_node(params_node, "threshold"))
22!
796
      threshold_ = std::stod(get_node_value(params_node, "threshold"));
22✔
797
    if (check_for_node(params_node, "ratio")) {
22!
798
      ratio_ = std::stod(get_node_value(params_node, "ratio"));
22✔
799
    }
800
  }
801

802
  // check update parameter values
803
  if (tally_value_ != "mean" && tally_value_ != "rel_err") {
82!
804
    fatal_error(fmt::format("Unsupported tally value '{}' specified for "
×
805
                            "weight window generation.",
806
      tally_value_));
×
807
  }
808
  if (threshold_ <= 0.0)
82!
809
    fatal_error(fmt::format("Invalid relative error threshold '{}' (<= 0.0) "
×
810
                            "specified for weight window generation",
811
      ratio_));
×
812
  if (ratio_ <= 1.0)
82!
813
    fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) "
×
814
                            "specified for weight window generation"));
815

816
  // create a matching weight windows object
817
  auto wws = WeightWindows::create();
82✔
818
  ww_idx_ = wws->index();
82✔
819
  wws->set_mesh(mesh_idx);
82✔
820
  if (e_bounds.size() > 0)
82!
821
    wws->set_energy_bounds(e_bounds);
82✔
822
  wws->set_particle_type(particle_type);
82✔
823
  wws->set_defaults();
82✔
824
}
82✔
825

826
void WeightWindowsGenerator::create_tally()
82✔
827
{
828
  const auto& wws = variance_reduction::weight_windows[ww_idx_];
82✔
829

830
  // create a tally based on the WWG information
831
  Tally* ww_tally = Tally::create();
82✔
832
  tally_idx_ = model::tally_map[ww_tally->id()];
82✔
833
  ww_tally->set_scores({"flux"});
164!
834

835
  int32_t mesh_id = wws->mesh()->id();
82✔
836
  int32_t mesh_idx = model::mesh_map.at(mesh_id);
82✔
837
  // see if there's already a mesh filter using this mesh
838
  bool found_mesh_filter = false;
82✔
839
  for (const auto& f : model::tally_filters) {
259✔
840
    if (f->type() == FilterType::MESH) {
188✔
841
      const auto* mesh_filter = dynamic_cast<MeshFilter*>(f.get());
11!
842
      if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated() &&
22!
843
          !mesh_filter->rotated()) {
11!
844
        ww_tally->add_filter(f.get());
11✔
845
        found_mesh_filter = true;
11✔
846
        break;
11✔
847
      }
848
    }
849
  }
850

851
  if (!found_mesh_filter) {
82✔
852
    auto mesh_filter = Filter::create("mesh");
71✔
853
    openmc_mesh_filter_set_mesh(mesh_filter->index(), model::mesh_map[mesh_id]);
71✔
854
    ww_tally->add_filter(mesh_filter);
71✔
855
  }
856

857
  const auto& e_bounds = wws->energy_bounds();
82✔
858
  if (e_bounds.size() > 0) {
82!
859
    auto energy_filter = Filter::create("energy");
82✔
860
    openmc_energy_filter_set_bins(
164✔
861
      energy_filter->index(), e_bounds.size(), e_bounds.data());
82✔
862
    ww_tally->add_filter(energy_filter);
82✔
863
  }
864

865
  // add a particle filter
866
  auto particle_type = wws->particle_type();
82✔
867
  auto particle_filter = Filter::create("particle");
82✔
868
  auto pf = dynamic_cast<ParticleFilter*>(particle_filter);
82!
869
  pf->set_particles({&particle_type, 1});
82✔
870
  ww_tally->add_filter(particle_filter);
82✔
871
}
82✔
872

873
void WeightWindowsGenerator::update() const
2,417✔
874
{
875
  const auto& wws = variance_reduction::weight_windows[ww_idx_];
2,417✔
876

877
  Tally* tally = model::tallies[tally_idx_].get();
2,417✔
878

879
  // If in random ray mode, only update on the last batch
880
  if (settings::solver_type == SolverType::RANDOM_RAY) {
2,417✔
881
    if (simulation::current_batch != settings::n_batches) {
2,252✔
882
      return;
2,154✔
883
    }
884
    // If in Monte Carlo mode and beyond the number of max realizations or
885
    // not at the correct update interval, skip the update
886
  } else if (max_realizations_ < tally->n_realizations_ ||
165✔
887
             tally->n_realizations_ % update_interval_ != 0) {
33!
888
    return;
132✔
889
  }
890

891
  wws->update_weights(tally, tally_value_, threshold_, ratio_, method_);
131✔
892

893
  // if we're not doing on the fly generation, reset the tally results once
894
  // we're done with the update
895
  if (!on_the_fly_)
131!
896
    tally->reset();
×
897

898
  // TODO: deactivate or remove tally once weight window generation is
899
  // complete
900
}
901

902
//==============================================================================
903
// Non-member functions
904
//==============================================================================
905

906
WeightWindow search_weight_window(const Particle& p)
82,252,639✔
907
{
908
  // TODO: this is a linear search - should do something more clever
909
  int mesh_bin;
910
  WeightWindow weight_window;
82,252,639✔
911
  for (const auto& ww : variance_reduction::weight_windows) {
115,666,136✔
912
    mesh_bin = ww->get_mesh_bin(p);
90,902,269✔
913
    if (mesh_bin < 0)
90,902,269✔
914
      continue;
19,817,025✔
915
    weight_window = ww->get_weight_window(p.E(), mesh_bin);
71,085,244✔
916
    if (weight_window.is_valid())
71,085,244✔
917
      return weight_window;
57,488,772✔
918
  }
919
  return {};
24,763,867✔
920
}
921

922
void apply_weight_windows(Particle& p)
167,087,795✔
923
{
924
  if (!settings::weight_windows_on)
167,087,795✔
925
    return;
166,343,054✔
926

927
  // WW on photon and neutron only
928
  if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon)
748,418!
NEW
929
    return;
×
930

931
  // skip dead or no energy
932
  if (p.E() <= 0 || !p.alive())
748,418!
933
    return;
3,677✔
934

935
  auto ww = search_weight_window(p);
744,741✔
936
  if (ww.is_valid()) {
744,741✔
937
    apply_weight_window(p, ww);
677,508✔
938
  } else {
939
    if (p.wgt_ww_born() == -1.0)
67,233✔
940
      p.wgt_ww_born() = 1.0;
65,802✔
941
  }
942
}
943

944
void apply_weight_window(Particle& p, WeightWindow weight_window)
57,488,772✔
945
{
946
  // skip dead or no energy
947
  if (p.E() <= 0 || !p.alive())
57,488,772✔
948
    return;
3,596,046✔
949

950
  // If particle has not yet had its birth weight window value set, set it to
951
  // the current weight window.
952
  if (p.wgt_ww_born() == -1.0)
53,892,726✔
953
    p.wgt_ww_born() =
671,878✔
954
      (weight_window.lower_weight + weight_window.upper_weight) / 2;
671,878✔
955

956
  // Normalize weight windows based on particle's starting weight
957
  // and the value of the weight window the particle was born in.
958
  weight_window.scale(p.wgt_born() / p.wgt_ww_born());
53,892,726✔
959

960
  // get the paramters
961
  double weight = p.wgt();
53,892,726✔
962

963
  // first check to see if particle should be killed for weight cutoff
964
  if (p.wgt() < weight_window.weight_cutoff) {
53,892,726!
NEW
965
    p.wgt() = 0.0;
×
NEW
966
    return;
×
967
  }
968

969
  // check if particle is far above current weight window
970
  // only do this if the factor is not already set on the particle and a
971
  // maximum lower bound ratio is specified
972
  if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 &&
53,895,564✔
973
      p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) {
2,838!
974
    p.ww_factor() =
2,838✔
975
      p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio);
2,838✔
976
  }
977

978
  // move weight window closer to the particle weight if needed
979
  if (p.ww_factor() > 1.0)
53,892,726✔
980
    weight_window.scale(p.ww_factor());
1,356,443✔
981

982
  // if particle's weight is above the weight window split until they are within
983
  // the window
984
  if (weight > weight_window.upper_weight) {
53,892,726✔
985
    // do not further split the particle if above the limit
986
    if (p.n_split() >= settings::max_history_splits)
13,452,126✔
987
      return;
12,139,997✔
988

989
    double n_split = std::ceil(weight / weight_window.upper_weight);
1,312,129✔
990
    double max_split = weight_window.max_split;
1,312,129✔
991
    n_split = std::min(n_split, max_split);
1,312,129✔
992

993
    p.n_split() += n_split;
1,312,129✔
994

995
    // Create secondaries and divide weight among all particles
996
    int i_split = std::round(n_split);
1,312,129✔
997
    for (int l = 0; l < i_split - 1; l++) {
5,399,706✔
998
      p.split(weight / n_split);
4,087,577✔
999
    }
1000
    // remaining weight is applied to current particle
1001
    p.wgt() = weight / n_split;
1,312,129✔
1002

1003
  } else if (weight <= weight_window.lower_weight) {
40,440,600✔
1004
    // if the particle weight is below the window, play Russian roulette
1005
    double weight_survive =
1006
      std::min(weight * weight_window.max_split, weight_window.survival_weight);
1,268,057✔
1007
    russian_roulette(p, weight_survive);
1,268,057✔
1008
  } // else particle is in the window, continue as normal
1009
}
1010

1011
void free_memory_weight_windows()
8,153✔
1012
{
1013
  variance_reduction::ww_map.clear();
8,153✔
1014
  variance_reduction::weight_windows.clear();
8,153✔
1015
}
8,153✔
1016

1017
void finalize_variance_reduction()
8,006✔
1018
{
1019
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
8,088✔
1020
    wwg->create_tally();
82✔
1021
  }
1022
}
8,006✔
1023

1024
//==============================================================================
1025
// C API
1026
//==============================================================================
1027

1028
int verify_ww_index(int32_t index)
1,991✔
1029
{
1030
  if (index < 0 || index >= variance_reduction::weight_windows.size()) {
1,991!
1031
    set_errmsg(fmt::format("Index '{}' for weight windows is invalid", index));
×
1032
    return OPENMC_E_OUT_OF_BOUNDS;
×
1033
  }
1034
  return 0;
1,991✔
1035
}
1036

1037
extern "C" int openmc_get_weight_windows_index(int32_t id, int32_t* idx)
165✔
1038
{
1039
  auto it = variance_reduction::ww_map.find(id);
165✔
1040
  if (it == variance_reduction::ww_map.end()) {
165!
1041
    set_errmsg(fmt::format("No weight windows exist with ID={}", id));
×
1042
    return OPENMC_E_INVALID_ID;
×
1043
  }
1044

1045
  *idx = it->second;
165✔
1046
  return 0;
165✔
1047
}
1048

1049
extern "C" int openmc_weight_windows_get_id(int32_t index, int32_t* id)
517✔
1050
{
1051
  if (int err = verify_ww_index(index))
517!
1052
    return err;
×
1053

1054
  const auto& wws = variance_reduction::weight_windows.at(index);
517✔
1055
  *id = wws->id();
517✔
1056
  return 0;
517✔
1057
}
1058

1059
extern "C" int openmc_weight_windows_set_id(int32_t index, int32_t id)
154✔
1060
{
1061
  if (int err = verify_ww_index(index))
154!
1062
    return err;
×
1063

1064
  const auto& wws = variance_reduction::weight_windows.at(index);
154✔
1065
  wws->set_id(id);
154✔
1066
  return 0;
154✔
1067
}
1068

1069
extern "C" int openmc_weight_windows_update_magic(int32_t ww_idx,
121✔
1070
  int32_t tally_idx, const char* value, double threshold, double ratio)
1071
{
1072
  if (int err = verify_ww_index(ww_idx))
121!
1073
    return err;
×
1074

1075
  if (tally_idx < 0 || tally_idx >= model::tallies.size()) {
121!
1076
    set_errmsg(fmt::format("Index '{}' for tally is invalid", tally_idx));
×
1077
    return OPENMC_E_OUT_OF_BOUNDS;
×
1078
  }
1079

1080
  // get the requested tally
1081
  const Tally* tally = model::tallies.at(tally_idx).get();
121✔
1082

1083
  // get the WeightWindows object
1084
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
121✔
1085

1086
  wws->update_weights(tally, value, threshold, ratio);
121✔
1087

1088
  return 0;
121✔
1089
}
1090

1091
extern "C" int openmc_weight_windows_set_mesh(int32_t ww_idx, int32_t mesh_idx)
154✔
1092
{
1093
  if (int err = verify_ww_index(ww_idx))
154!
1094
    return err;
×
1095
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
154✔
1096
  wws->set_mesh(mesh_idx);
154✔
1097
  return 0;
154✔
1098
}
1099

1100
extern "C" int openmc_weight_windows_get_mesh(int32_t ww_idx, int32_t* mesh_idx)
11✔
1101
{
1102
  if (int err = verify_ww_index(ww_idx))
11!
1103
    return err;
×
1104
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
11✔
1105
  *mesh_idx = model::mesh_map.at(wws->mesh()->id());
11✔
1106
  return 0;
11✔
1107
}
1108

1109
extern "C" int openmc_weight_windows_set_energy_bounds(
132✔
1110
  int32_t ww_idx, double* e_bounds, size_t e_bounds_size)
1111
{
1112
  if (int err = verify_ww_index(ww_idx))
132!
1113
    return err;
×
1114
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
132✔
1115
  wws->set_energy_bounds({e_bounds, e_bounds_size});
132✔
1116
  return 0;
132✔
1117
}
1118

1119
extern "C" int openmc_weight_windows_get_energy_bounds(
11✔
1120
  int32_t ww_idx, const double** e_bounds, size_t* e_bounds_size)
1121
{
1122
  if (int err = verify_ww_index(ww_idx))
11!
1123
    return err;
×
1124
  const auto& wws = variance_reduction::weight_windows[ww_idx].get();
11✔
1125
  *e_bounds = wws->energy_bounds().data();
11✔
1126
  *e_bounds_size = wws->energy_bounds().size();
11✔
1127
  return 0;
11✔
1128
}
1129

1130
extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle)
176✔
1131
{
1132
  if (int err = verify_ww_index(index))
176!
1133
    return err;
×
1134

1135
  const auto& wws = variance_reduction::weight_windows.at(index);
176✔
1136
  wws->set_particle_type(static_cast<ParticleType>(particle));
176✔
1137
  return 0;
176✔
1138
}
1139

1140
extern "C" int openmc_weight_windows_get_particle(int32_t index, int* particle)
44✔
1141
{
1142
  if (int err = verify_ww_index(index))
44!
1143
    return err;
×
1144

1145
  const auto& wws = variance_reduction::weight_windows.at(index);
44✔
1146
  *particle = static_cast<int>(wws->particle_type());
44✔
1147
  return 0;
44✔
1148
}
1149

1150
extern "C" int openmc_weight_windows_get_bounds(int32_t index,
484✔
1151
  const double** lower_bounds, const double** upper_bounds, size_t* size)
1152
{
1153
  if (int err = verify_ww_index(index))
484!
1154
    return err;
×
1155

1156
  const auto& wws = variance_reduction::weight_windows[index];
484✔
1157
  *size = wws->lower_ww_bounds().size();
484✔
1158
  *lower_bounds = wws->lower_ww_bounds().data();
484✔
1159
  *upper_bounds = wws->upper_ww_bounds().data();
484✔
1160
  return 0;
484✔
1161
}
1162

1163
extern "C" int openmc_weight_windows_set_bounds(int32_t index,
11✔
1164
  const double* lower_bounds, const double* upper_bounds, size_t size)
1165
{
1166
  if (int err = verify_ww_index(index))
11!
1167
    return err;
×
1168

1169
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1170
  wws->set_bounds({lower_bounds, size}, {upper_bounds, size});
11✔
1171
  return 0;
11✔
1172
}
1173

1174
extern "C" int openmc_weight_windows_get_survival_ratio(
33✔
1175
  int32_t index, double* ratio)
1176
{
1177
  if (int err = verify_ww_index(index))
33!
1178
    return err;
×
1179
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1180
  *ratio = wws->survival_ratio();
33✔
1181
  return 0;
33✔
1182
}
1183

1184
extern "C" int openmc_weight_windows_set_survival_ratio(
11✔
1185
  int32_t index, double ratio)
1186
{
1187
  if (int err = verify_ww_index(index))
11!
1188
    return err;
×
1189
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1190
  wws->survival_ratio() = ratio;
11✔
1191
  std::cout << "Survival ratio: " << wws->survival_ratio() << std::endl;
11✔
1192
  return 0;
11✔
1193
}
1194

1195
extern "C" int openmc_weight_windows_get_max_lower_bound_ratio(
33✔
1196
  int32_t index, double* lb_ratio)
1197
{
1198
  if (int err = verify_ww_index(index))
33!
1199
    return err;
×
1200
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1201
  *lb_ratio = wws->max_lower_bound_ratio();
33✔
1202
  return 0;
33✔
1203
}
1204

1205
extern "C" int openmc_weight_windows_set_max_lower_bound_ratio(
11✔
1206
  int32_t index, double lb_ratio)
1207
{
1208
  if (int err = verify_ww_index(index))
11!
1209
    return err;
×
1210
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1211
  wws->max_lower_bound_ratio() = lb_ratio;
11✔
1212
  return 0;
11✔
1213
}
1214

1215
extern "C" int openmc_weight_windows_get_weight_cutoff(
33✔
1216
  int32_t index, double* cutoff)
1217
{
1218
  if (int err = verify_ww_index(index))
33!
1219
    return err;
×
1220
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1221
  *cutoff = wws->weight_cutoff();
33✔
1222
  return 0;
33✔
1223
}
1224

1225
extern "C" int openmc_weight_windows_set_weight_cutoff(
11✔
1226
  int32_t index, double cutoff)
1227
{
1228
  if (int err = verify_ww_index(index))
11!
1229
    return err;
×
1230
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1231
  wws->weight_cutoff() = cutoff;
11✔
1232
  return 0;
11✔
1233
}
1234

1235
extern "C" int openmc_weight_windows_get_max_split(
33✔
1236
  int32_t index, int* max_split)
1237
{
1238
  if (int err = verify_ww_index(index))
33!
1239
    return err;
×
1240
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1241
  *max_split = wws->max_split();
33✔
1242
  return 0;
33✔
1243
}
1244

1245
extern "C" int openmc_weight_windows_set_max_split(int32_t index, int max_split)
11✔
1246
{
1247
  if (int err = verify_ww_index(index))
11!
1248
    return err;
×
1249
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1250
  wws->max_split() = max_split;
11✔
1251
  return 0;
11✔
1252
}
1253

1254
extern "C" int openmc_extend_weight_windows(
154✔
1255
  int32_t n, int32_t* index_start, int32_t* index_end)
1256
{
1257
  if (index_start)
154!
1258
    *index_start = variance_reduction::weight_windows.size();
154✔
1259
  if (index_end)
154!
1260
    *index_end = variance_reduction::weight_windows.size() + n - 1;
×
1261
  for (int i = 0; i < n; ++i)
308✔
1262
    variance_reduction::weight_windows.push_back(make_unique<WeightWindows>());
154✔
1263
  return 0;
154✔
1264
}
1265

1266
extern "C" size_t openmc_weight_windows_size()
154✔
1267
{
1268
  return variance_reduction::weight_windows.size();
154✔
1269
}
1270

1271
extern "C" int openmc_weight_windows_export(const char* filename)
164✔
1272
{
1273

1274
  if (!mpi::master)
164✔
1275
    return 0;
30✔
1276

1277
  std::string name = filename ? filename : "weight_windows.h5";
268✔
1278

1279
  write_message(fmt::format("Exporting weight windows to {}...", name), 5);
134✔
1280

1281
  hid_t ww_file = file_open(name, 'w');
134✔
1282

1283
  // Write file type
1284
  write_attribute(ww_file, "filetype", "weight_windows");
134✔
1285

1286
  // Write revisiion number for state point file
1287
  write_attribute(ww_file, "version", VERSION_WEIGHT_WINDOWS);
134✔
1288

1289
  hid_t weight_windows_group = create_group(ww_file, "weight_windows");
134✔
1290

1291
  hid_t mesh_group = create_group(ww_file, "meshes");
134✔
1292

1293
  std::vector<int32_t> mesh_ids;
134✔
1294
  std::vector<int32_t> ww_ids;
134✔
1295
  for (const auto& ww : variance_reduction::weight_windows) {
268✔
1296

1297
    ww->to_hdf5(weight_windows_group);
134✔
1298
    ww_ids.push_back(ww->id());
134✔
1299

1300
    // if the mesh has already been written, move on
1301
    int32_t mesh_id = ww->mesh()->id();
134✔
1302
    if (std::find(mesh_ids.begin(), mesh_ids.end(), mesh_id) != mesh_ids.end())
134!
1303
      continue;
×
1304

1305
    mesh_ids.push_back(mesh_id);
134✔
1306
    ww->mesh()->to_hdf5(mesh_group);
134✔
1307
  }
1308

1309
  write_attribute(mesh_group, "n_meshes", mesh_ids.size());
134✔
1310
  write_attribute(mesh_group, "ids", mesh_ids);
134✔
1311
  close_group(mesh_group);
134✔
1312

1313
  write_attribute(weight_windows_group, "n_weight_windows", ww_ids.size());
134✔
1314
  write_attribute(weight_windows_group, "ids", ww_ids);
134✔
1315
  close_group(weight_windows_group);
134✔
1316

1317
  file_close(ww_file);
134✔
1318

1319
  return 0;
134✔
1320
}
134✔
1321

1322
extern "C" int openmc_weight_windows_import(const char* filename)
11✔
1323
{
1324
  std::string name = filename ? filename : "weight_windows.h5";
11!
1325

1326
  if (mpi::master)
11!
1327
    write_message(fmt::format("Importing weight windows from {}...", name), 5);
11✔
1328

1329
  if (!file_exists(name)) {
11!
1330
    set_errmsg(fmt::format("File '{}' does not exist", name));
×
1331
  }
1332

1333
  hid_t ww_file = file_open(name, 'r');
11✔
1334

1335
  // Check that filetype is correct
1336
  std::string filetype;
11✔
1337
  read_attribute(ww_file, "filetype", filetype);
11✔
1338
  if (filetype != "weight_windows") {
11!
1339
    file_close(ww_file);
×
1340
    set_errmsg(fmt::format("File '{}' is not a weight windows file.", name));
×
1341
    return OPENMC_E_INVALID_ARGUMENT;
×
1342
  }
1343

1344
  // Check that the file version is compatible
1345
  std::array<int, 2> file_version;
1346
  read_attribute(ww_file, "version", file_version);
11✔
1347
  if (file_version[0] != VERSION_WEIGHT_WINDOWS[0]) {
11!
1348
    std::string err_msg =
1349
      fmt::format("File '{}' has version {} which is incompatible with the "
1350
                  "expected version ({}).",
1351
        name, file_version, VERSION_WEIGHT_WINDOWS);
×
1352
    set_errmsg(err_msg);
×
1353
    return OPENMC_E_INVALID_ARGUMENT;
×
1354
  }
×
1355

1356
  hid_t weight_windows_group = open_group(ww_file, "weight_windows");
11✔
1357

1358
  hid_t mesh_group = open_group(ww_file, "meshes");
11✔
1359

1360
  read_meshes(mesh_group);
11✔
1361

1362
  std::vector<std::string> names = group_names(weight_windows_group);
11✔
1363

1364
  for (const auto& name : names) {
22✔
1365
    WeightWindows::from_hdf5(weight_windows_group, name);
11✔
1366
  }
1367

1368
  close_group(weight_windows_group);
11✔
1369

1370
  file_close(ww_file);
11✔
1371

1372
  return 0;
11✔
1373
}
11✔
1374

1375
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc