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

openmc-dev / openmc / 20472454196

23 Dec 2025 09:56PM UTC coverage: 82.155% (+0.02%) from 82.139%
20472454196

Pull #3692

github

web-flow
Merge 064b0c3b2 into 3f06a42ab
Pull Request #3692: Fix a bug in rotational periodic boundary conditions

17067 of 23647 branches covered (72.17%)

Branch coverage included in aggregate %.

29 of 33 new or added lines in 3 files covered. (87.88%)

185 existing lines in 9 files now uncovered.

55183 of 64296 relevant lines covered (85.83%)

43492701.44 hits per line

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

78.87
/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
// Non-member functions
54
//==============================================================================
55

56
void apply_weight_windows(Particle& p)
2,147,483,647✔
57
{
58
  if (!settings::weight_windows_on)
2,147,483,647✔
59
    return;
2,147,483,647✔
60

61
  // WW on photon and neutron only
62
  if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon)
82,243,829✔
63
    return;
10,928,406✔
64

65
  // skip dead or no energy
66
  if (p.E() <= 0 || !p.alive())
71,315,423✔
67
    return;
4,092,414✔
68

69
  bool in_domain = false;
67,223,009✔
70
  // TODO: this is a linear search - should do something more clever
71
  WeightWindow weight_window;
67,223,009✔
72
  for (const auto& ww : variance_reduction::weight_windows) {
84,612,771✔
73
    weight_window = ww->get_weight_window(p);
71,271,625✔
74
    if (weight_window.is_valid())
71,271,625✔
75
      break;
53,881,863✔
76
  }
77

78
  // If particle has not yet had its birth weight window value set, set it to
79
  // the current weight window (or 1.0 if not born in a weight window).
80
  if (p.wgt_ww_born() == -1.0) {
67,223,009✔
81
    if (weight_window.is_valid()) {
737,680✔
82
      p.wgt_ww_born() =
671,878✔
83
        (weight_window.lower_weight + weight_window.upper_weight) / 2;
671,878✔
84
    } else {
85
      p.wgt_ww_born() = 1.0;
65,802✔
86
    }
87
  }
88

89
  // particle is not in any of the ww domains, do nothing
90
  if (!weight_window.is_valid())
67,223,009✔
91
    return;
13,341,146✔
92

93
  // Normalize weight windows based on particle's starting weight
94
  // and the value of the weight window the particle was born in.
95
  weight_window.scale(p.wgt_born() / p.wgt_ww_born());
53,881,863✔
96

97
  // get the paramters
98
  double weight = p.wgt();
53,881,863✔
99

100
  // first check to see if particle should be killed for weight cutoff
101
  if (p.wgt() < weight_window.weight_cutoff) {
53,881,863!
102
    p.wgt() = 0.0;
×
103
    return;
×
104
  }
105

106
  // check if particle is far above current weight window
107
  // only do this if the factor is not already set on the particle and a
108
  // maximum lower bound ratio is specified
109
  if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 &&
53,884,701✔
110
      p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) {
2,838!
111
    p.ww_factor() =
2,838✔
112
      p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio);
2,838✔
113
  }
114

115
  // move weight window closer to the particle weight if needed
116
  if (p.ww_factor() > 1.0)
53,881,863✔
117
    weight_window.scale(p.ww_factor());
1,356,443✔
118

119
  // if particle's weight is above the weight window split until they are within
120
  // the window
121
  if (weight > weight_window.upper_weight) {
53,881,863✔
122
    // do not further split the particle if above the limit
123
    if (p.n_split() >= settings::max_history_splits)
13,447,769✔
124
      return;
12,135,589✔
125

126
    double n_split = std::ceil(weight / weight_window.upper_weight);
1,312,180✔
127
    double max_split = weight_window.max_split;
1,312,180✔
128
    n_split = std::min(n_split, max_split);
1,312,180✔
129

130
    p.n_split() += n_split;
1,312,180✔
131

132
    // Create secondaries and divide weight among all particles
133
    int i_split = std::round(n_split);
1,312,180✔
134
    for (int l = 0; l < i_split - 1; l++) {
5,400,118✔
135
      p.split(weight / n_split);
4,087,938✔
136
    }
137
    // remaining weight is applied to current particle
138
    p.wgt() = weight / n_split;
1,312,180✔
139

140
  } else if (weight <= weight_window.lower_weight) {
40,434,094✔
141
    // if the particle weight is below the window, play Russian roulette
142
    double weight_survive =
143
      std::min(weight * weight_window.max_split, weight_window.survival_weight);
1,268,044✔
144
    russian_roulette(p, weight_survive);
1,268,044✔
145
  } // else particle is in the window, continue as normal
146
}
147

148
void free_memory_weight_windows()
7,771✔
149
{
150
  variance_reduction::ww_map.clear();
7,771✔
151
  variance_reduction::weight_windows.clear();
7,771✔
152
}
7,771✔
153

154
//==============================================================================
155
// WeightWindowSettings implementation
156
//==============================================================================
157

158
WeightWindows::WeightWindows(int32_t id)
247✔
159
{
160
  index_ = variance_reduction::weight_windows.size();
247✔
161
  set_id(id);
247✔
162
  set_defaults();
247✔
163
}
247✔
164

165
WeightWindows::WeightWindows(pugi::xml_node node)
92✔
166
{
167
  // Make sure required elements are present
168
  const vector<std::string> required_elems {
169
    "id", "particle_type", "lower_ww_bounds", "upper_ww_bounds"};
644✔
170
  for (const auto& elem : required_elems) {
460✔
171
    if (!check_for_node(node, elem.c_str())) {
368!
172
      fatal_error(fmt::format("Must specify <{}> for weight windows.", elem));
×
173
    }
174
  }
175

176
  // Get weight windows ID
177
  int32_t id = std::stoi(get_node_value(node, "id"));
92✔
178
  this->set_id(id);
92✔
179

180
  // get the particle type
181
  auto particle_type_str = std::string(get_node_value(node, "particle_type"));
92✔
182
  particle_type_ = openmc::str_to_particle_type(particle_type_str);
92✔
183

184
  // Determine associated mesh
185
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
92✔
186
  set_mesh(model::mesh_map.at(mesh_id));
92✔
187

188
  // energy bounds
189
  if (check_for_node(node, "energy_bounds"))
92✔
190
    energy_bounds_ = get_node_array<double>(node, "energy_bounds");
77✔
191

192
  // get the survival value - optional
193
  if (check_for_node(node, "survival_ratio")) {
92!
194
    survival_ratio_ = std::stod(get_node_value(node, "survival_ratio"));
92✔
195
    if (survival_ratio_ <= 1)
92!
196
      fatal_error("Survival to lower weight window ratio must bigger than 1 "
×
197
                  "and less than the upper to lower weight window ratio.");
198
  }
199

200
  // get the max lower bound ratio - optional
201
  if (check_for_node(node, "max_lower_bound_ratio")) {
92✔
202
    max_lb_ratio_ = std::stod(get_node_value(node, "max_lower_bound_ratio"));
34✔
203
    if (max_lb_ratio_ < 1.0) {
34!
204
      fatal_error("Maximum lower bound ratio must be larger than 1");
×
205
    }
206
  }
207

208
  // get the max split - optional
209
  if (check_for_node(node, "max_split")) {
92!
210
    max_split_ = std::stod(get_node_value(node, "max_split"));
92✔
211
    if (max_split_ <= 1)
92!
212
      fatal_error("max split must be larger than 1");
×
213
  }
214

215
  // weight cutoff - optional
216
  if (check_for_node(node, "weight_cutoff")) {
92!
217
    weight_cutoff_ = std::stod(get_node_value(node, "weight_cutoff"));
92✔
218
    if (weight_cutoff_ <= 0)
92!
219
      fatal_error("weight_cutoff must be larger than 0");
×
220
    if (weight_cutoff_ > 1)
92!
221
      fatal_error("weight_cutoff must be less than 1");
×
222
  }
223

224
  // read the lower/upper weight bounds
225
  this->set_bounds(get_node_array<double>(node, "lower_ww_bounds"),
92✔
226
    get_node_array<double>(node, "upper_ww_bounds"));
184✔
227

228
  set_defaults();
92✔
229
}
92✔
230

231
WeightWindows::~WeightWindows()
339✔
232
{
233
  variance_reduction::ww_map.erase(id());
339✔
234
}
339✔
235

236
WeightWindows* WeightWindows::create(int32_t id)
93✔
237
{
238
  variance_reduction::weight_windows.push_back(make_unique<WeightWindows>());
93✔
239
  auto wws = variance_reduction::weight_windows.back().get();
93✔
240
  variance_reduction::ww_map[wws->id()] =
93✔
241
    variance_reduction::weight_windows.size() - 1;
93✔
242
  return wws;
93✔
243
}
244

245
WeightWindows* WeightWindows::from_hdf5(
11✔
246
  hid_t wws_group, const std::string& group_name)
247
{
248
  // collect ID from the name of this group
249
  hid_t ww_group = open_group(wws_group, group_name);
11✔
250

251
  auto wws = WeightWindows::create();
11✔
252

253
  std::string particle_type;
11✔
254
  read_dataset(ww_group, "particle_type", particle_type);
11✔
255
  wws->particle_type_ = openmc::str_to_particle_type(particle_type);
11✔
256

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

259
  int32_t mesh_id;
260
  read_dataset(ww_group, "mesh", mesh_id);
11✔
261

262
  if (model::mesh_map.count(mesh_id) == 0) {
11!
263
    fatal_error(
×
264
      fmt::format("Mesh {} used in weight windows does not exist.", mesh_id));
×
265
  }
266
  wws->set_mesh(model::mesh_map[mesh_id]);
11✔
267

268
  wws->lower_ww_ = xt::empty<double>(wws->bounds_size());
11✔
269
  wws->upper_ww_ = xt::empty<double>(wws->bounds_size());
11✔
270

271
  read_dataset<double>(ww_group, "lower_ww_bounds", wws->lower_ww_);
11✔
272
  read_dataset<double>(ww_group, "upper_ww_bounds", wws->upper_ww_);
11✔
273
  read_dataset(ww_group, "survival_ratio", wws->survival_ratio_);
11✔
274
  read_dataset(ww_group, "max_lower_bound_ratio", wws->max_lb_ratio_);
11✔
275
  read_dataset(ww_group, "max_split", wws->max_split_);
11✔
276
  read_dataset(ww_group, "weight_cutoff", wws->weight_cutoff_);
11✔
277

278
  close_group(ww_group);
11✔
279

280
  return wws;
11✔
281
}
11✔
282

283
void WeightWindows::set_defaults()
421✔
284
{
285
  // set energy bounds to the min/max energy supported by the data
286
  if (energy_bounds_.size() == 0) {
421✔
287
    int p_type = static_cast<int>(particle_type_);
262✔
288
    energy_bounds_.push_back(data::energy_min[p_type]);
262✔
289
    energy_bounds_.push_back(data::energy_max[p_type]);
262✔
290
  }
291
}
421✔
292

293
void WeightWindows::allocate_ww_bounds()
553✔
294
{
295
  auto shape = bounds_size();
553✔
296
  if (shape[0] * shape[1] == 0) {
553!
297
    auto msg = fmt::format(
298
      "Size of weight window bounds is zero for WeightWindows {}", id());
×
299
    warning(msg);
×
300
  }
×
301
  lower_ww_ = xt::empty<double>(shape);
553✔
302
  lower_ww_.fill(-1);
553✔
303
  upper_ww_ = xt::empty<double>(shape);
553✔
304
  upper_ww_.fill(-1);
553✔
305
}
553✔
306

307
void WeightWindows::set_id(int32_t id)
493✔
308
{
309
  assert(id >= 0 || id == C_NONE);
402!
310

311
  // Clear entry in mesh map in case one was already assigned
312
  if (id_ != C_NONE) {
493!
313
    variance_reduction::ww_map.erase(id_);
493✔
314
    id_ = C_NONE;
493✔
315
  }
316

317
  // Ensure no other mesh has the same ID
318
  if (variance_reduction::ww_map.find(id) != variance_reduction::ww_map.end()) {
493!
319
    throw std::runtime_error {
×
320
      fmt::format("Two weight windows have the same ID: {}", id)};
×
321
  }
322

323
  // If no ID is specified, auto-assign the next ID in the sequence
324
  if (id == C_NONE) {
493✔
325
    id = 0;
247✔
326
    for (const auto& m : variance_reduction::weight_windows) {
269✔
327
      id = std::max(id, m->id_);
22✔
328
    }
329
    ++id;
247✔
330
  }
331

332
  // Update ID and entry in the mesh map
333
  id_ = id;
493✔
334
  variance_reduction::ww_map[id] = index_;
493✔
335
}
493✔
336

337
void WeightWindows::set_energy_bounds(span<const double> bounds)
214✔
338
{
339
  energy_bounds_.clear();
214✔
340
  energy_bounds_.insert(energy_bounds_.begin(), bounds.begin(), bounds.end());
214✔
341
  // if the mesh is set, allocate space for weight window bounds
342
  if (mesh_idx_ != C_NONE)
214!
343
    allocate_ww_bounds();
214✔
344
}
214✔
345

346
void WeightWindows::set_particle_type(ParticleType p_type)
258✔
347
{
348
  if (p_type != ParticleType::neutron && p_type != ParticleType::photon)
258!
349
    fatal_error(
×
350
      fmt::format("Particle type '{}' cannot be applied to weight windows.",
×
351
        particle_type_to_str(p_type)));
×
352
  particle_type_ = p_type;
258✔
353
}
258✔
354

355
void WeightWindows::set_mesh(int32_t mesh_idx)
339✔
356
{
357
  if (mesh_idx < 0 || mesh_idx >= model::meshes.size())
339!
358
    fatal_error(fmt::format("Could not find a mesh for index {}", mesh_idx));
×
359

360
  mesh_idx_ = mesh_idx;
339✔
361
  model::meshes[mesh_idx_]->prepare_for_point_location();
339✔
362
  allocate_ww_bounds();
339✔
363
}
339✔
364

365
void WeightWindows::set_mesh(const std::unique_ptr<Mesh>& mesh)
×
366
{
367
  set_mesh(mesh.get());
×
368
}
×
369

370
void WeightWindows::set_mesh(const Mesh* mesh)
×
371
{
372
  set_mesh(model::mesh_map[mesh->id_]);
×
373
}
×
374

375
WeightWindow WeightWindows::get_weight_window(const Particle& p) const
71,271,625✔
376
{
377
  // check for particle type
378
  if (particle_type_ != p.type()) {
71,271,625✔
379
    return {};
3,872,836✔
380
  }
381

382
  // Get mesh index for particle's position
383
  const auto& mesh = this->mesh();
67,398,789✔
384
  int mesh_bin = mesh->get_bin(p.r());
67,398,789✔
385

386
  // particle is outside the weight window mesh
387
  if (mesh_bin < 0)
67,398,789✔
388
    return {};
14,588✔
389

390
  // particle energy
391
  double E = p.E();
67,384,201✔
392

393
  // check to make sure energy is in range, expects sorted energy values
394
  if (E < energy_bounds_.front() || E > energy_bounds_.back())
67,384,201!
395
    return {};
92,598✔
396

397
  // get the mesh bin in energy group
398
  int energy_bin =
399
    lower_bound_index(energy_bounds_.begin(), energy_bounds_.end(), E);
67,291,603✔
400

401
  // mesh_bin += energy_bin * mesh->n_bins();
402
  // Create individual weight window
403
  WeightWindow ww;
67,291,603✔
404
  ww.lower_weight = lower_ww_(energy_bin, mesh_bin);
67,291,603✔
405
  ww.upper_weight = upper_ww_(energy_bin, mesh_bin);
67,291,603✔
406
  ww.survival_weight = ww.lower_weight * survival_ratio_;
67,291,603✔
407
  ww.max_lb_ratio = max_lb_ratio_;
67,291,603✔
408
  ww.max_split = max_split_;
67,291,603✔
409
  ww.weight_cutoff = weight_cutoff_;
67,291,603✔
410
  return ww;
67,291,603✔
411
}
412

413
std::array<int, 2> WeightWindows::bounds_size() const
781✔
414
{
415
  int num_spatial_bins = this->mesh()->n_bins();
781✔
416
  int num_energy_bins =
417
    energy_bounds_.size() > 0 ? energy_bounds_.size() - 1 : 1;
781✔
418
  return {num_energy_bins, num_spatial_bins};
781✔
419
}
420

421
template<class T>
422
void WeightWindows::check_bounds(const T& lower, const T& upper) const
103✔
423
{
424
  // make sure that the upper and lower bounds have the same size
425
  if (lower.size() != upper.size()) {
103!
426
    auto msg = fmt::format("The upper and lower weight window lengths do not "
×
427
                           "match.\n Lower size: {}\n Upper size: {}",
428
      lower.size(), upper.size());
×
429
    fatal_error(msg);
×
430
  }
×
431
  this->check_bounds(lower);
103✔
432
}
103✔
433

434
template<class T>
435
void WeightWindows::check_bounds(const T& bounds) const
103✔
436
{
437
  // check that the number of weight window entries is correct
438
  auto dims = this->bounds_size();
103!
439
  if (bounds.size() != dims[0] * dims[1]) {
103!
440
    auto err_msg =
×
441
      fmt::format("In weight window domain {} the number of spatial "
442
                  "energy/spatial bins ({}) does not match the number "
443
                  "of weight bins ({})",
444
        id_, dims, bounds.size());
×
445
    fatal_error(err_msg);
×
446
  }
×
447
}
103✔
448

449
void WeightWindows::set_bounds(const xt::xtensor<double, 2>& lower_bounds,
×
450
  const xt::xtensor<double, 2>& upper_bounds)
451
{
452

453
  this->check_bounds(lower_bounds, upper_bounds);
×
454

455
  // set new weight window values
456
  lower_ww_ = lower_bounds;
×
457
  upper_ww_ = upper_bounds;
×
458
}
×
459

460
void WeightWindows::set_bounds(
×
461
  const xt::xtensor<double, 2>& lower_bounds, double ratio)
462
{
463
  this->check_bounds(lower_bounds);
×
464

465
  // set new weight window values
466
  lower_ww_ = lower_bounds;
×
467
  upper_ww_ = lower_bounds;
×
468
  upper_ww_ *= ratio;
×
469
}
×
470

471
void WeightWindows::set_bounds(
103✔
472
  span<const double> lower_bounds, span<const double> upper_bounds)
473
{
474
  check_bounds(lower_bounds, upper_bounds);
103✔
475
  auto shape = this->bounds_size();
103✔
476
  lower_ww_ = xt::empty<double>(shape);
103✔
477
  upper_ww_ = xt::empty<double>(shape);
103✔
478

479
  // set new weight window values
480
  xt::view(lower_ww_, xt::all()) =
206✔
481
    xt::adapt(lower_bounds.data(), lower_ww_.shape());
309✔
482
  xt::view(upper_ww_, xt::all()) =
206✔
483
    xt::adapt(upper_bounds.data(), upper_ww_.shape());
309✔
484
}
103✔
485

486
void WeightWindows::set_bounds(span<const double> lower_bounds, double ratio)
×
487
{
488
  this->check_bounds(lower_bounds);
×
489

490
  auto shape = this->bounds_size();
×
491
  lower_ww_ = xt::empty<double>(shape);
×
492
  upper_ww_ = xt::empty<double>(shape);
×
493

494
  // set new weight window values
495
  xt::view(lower_ww_, xt::all()) =
×
496
    xt::adapt(lower_bounds.data(), lower_ww_.shape());
×
497
  xt::view(upper_ww_, xt::all()) =
×
498
    xt::adapt(lower_bounds.data(), upper_ww_.shape());
×
499
  upper_ww_ *= ratio;
×
500
}
×
501

502
void WeightWindows::update_weights(const Tally* tally, const std::string& value,
252✔
503
  double threshold, double ratio, WeightWindowUpdateMethod method)
504
{
505
  ///////////////////////////
506
  // Setup and checks
507
  ///////////////////////////
508
  this->check_tally_update_compatibility(tally);
252✔
509

510
  // Dimensions of weight window arrays
511
  int e_bins = lower_ww_.shape()[0];
252✔
512
  int64_t mesh_bins = lower_ww_.shape()[1];
252✔
513

514
  // Initialize weight window arrays to -1.0 by default
515
#pragma omp parallel for collapse(2) schedule(static)
140✔
516
  for (int e = 0; e < e_bins; e++) {
934✔
517
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
518
      lower_ww_(e, m) = -1.0;
1,189,249✔
519
      upper_ww_(e, m) = -1.0;
1,189,249✔
520
    }
521
  }
522

523
  // determine which value to use
524
  const std::set<std::string> allowed_values = {"mean", "rel_err"};
1,260✔
525
  if (allowed_values.count(value) == 0) {
252!
526
    fatal_error(fmt::format("Invalid value '{}' specified for weight window "
×
527
                            "generation. Must be one of: 'mean' or 'rel_err'",
528
      value));
529
  }
530

531
  // determine the index of the specified score
532
  int score_index = tally->score_index("flux");
252✔
533
  if (score_index == C_NONE) {
252!
534
    fatal_error(
×
535
      fmt::format("A 'flux' score required for weight window generation "
×
536
                  "is not present on tally {}.",
537
        tally->id()));
×
538
  }
539

540
  ///////////////////////////
541
  // Extract tally data
542
  //
543
  // At the end of this section, the mean and rel_err array
544
  // is a 2D view of tally data (n_e_groups, n_mesh_bins)
545
  //
546
  ///////////////////////////
547

548
  // build a shape for a view of the tally results, this will always be
549
  // dimension 5 (3 filter dimensions, 1 score dimension, 1 results dimension)
550
  // Look for the size of the last dimension of the results array
551
  const auto& results_arr = tally->results();
252✔
552
  const int results_dim = static_cast<int>(results_arr.shape()[2]);
252✔
553
  std::array<int, 5> shape = {1, 1, 1, tally->n_scores(), results_dim};
252✔
554

555
  // set the shape for the filters applied on the tally
556
  for (int i = 0; i < tally->filters().size(); i++) {
964✔
557
    const auto& filter = model::tally_filters[tally->filters(i)];
712✔
558
    shape[i] = filter->n_bins();
712✔
559
  }
560

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

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

567
  // assign other filter types to dummy positions if needed
568
  if (!tally->has_filter(FilterType::PARTICLE))
252✔
569
    filter_types.push_back(FilterType::PARTICLE);
22✔
570

571
  if (!tally->has_filter(FilterType::ENERGY))
252✔
572
    filter_types.push_back(FilterType::ENERGY);
22✔
573

574
  // particle axis mapping
575
  transpose[0] =
252✔
576
    std::find(filter_types.begin(), filter_types.end(), FilterType::PARTICLE) -
252✔
577
    filter_types.begin();
252✔
578

579
  // energy axis mapping
580
  transpose[1] =
252✔
581
    std::find(filter_types.begin(), filter_types.end(), FilterType::ENERGY) -
252✔
582
    filter_types.begin();
252✔
583

584
  // mesh axis mapping
585
  transpose[2] =
252✔
586
    std::find(filter_types.begin(), filter_types.end(), FilterType::MESH) -
252✔
587
    filter_types.begin();
252✔
588

589
  // get a fully reshaped view of the tally according to tally ordering of
590
  // filters
591
  auto tally_values = xt::reshape_view(results_arr, shape);
252✔
592

593
  // get a that is (particle, energy, mesh, scores, values)
594
  auto transposed_view = xt::transpose(tally_values, transpose);
252✔
595

596
  // determine the dimension and index of the particle data
597
  int particle_idx = 0;
252✔
598
  if (tally->has_filter(FilterType::PARTICLE)) {
252✔
599
    // get the particle filter
600
    auto pf = tally->get_filter<ParticleFilter>();
230✔
601
    const auto& particles = pf->particles();
230✔
602

603
    // find the index of the particle that matches these weight windows
604
    auto p_it =
605
      std::find(particles.begin(), particles.end(), this->particle_type_);
230✔
606
    // if the particle filter doesn't have particle data for the particle
607
    // used on this weight windows instance, report an error
608
    if (p_it == particles.end()) {
230!
609
      auto msg = fmt::format("Particle type '{}' not present on Filter {} for "
610
                             "Tally {} used to update WeightWindows {}",
611
        particle_type_to_str(this->particle_type_), pf->id(), tally->id(),
×
612
        this->id());
×
613
      fatal_error(msg);
×
614
    }
×
615

616
    // use the index of the particle in the filter to down-select data later
617
    particle_idx = p_it - particles.begin();
230✔
618
  }
619

620
  // down-select data based on particle and score
621
  auto sum = xt::dynamic_view(
1,260✔
622
    transposed_view, {particle_idx, xt::all(), xt::all(), score_index,
504✔
623
                       static_cast<int>(TallyResult::SUM)});
1,008✔
624
  auto sum_sq = xt::dynamic_view(
1,260✔
625
    transposed_view, {particle_idx, xt::all(), xt::all(), score_index,
504✔
626
                       static_cast<int>(TallyResult::SUM_SQ)});
1,008✔
627
  int n = tally->n_realizations_;
252✔
628

629
  //////////////////////////////////////////////
630
  //
631
  // Assign new weight windows
632
  //
633
  // Use references to the existing weight window data
634
  // to store and update the values
635
  //
636
  //////////////////////////////////////////////
637

638
  // up to this point the data arrays are views into the tally results (no
639
  // computation has been performed) now we'll switch references to the tally's
640
  // bounds to avoid allocating additional memory
641
  auto& new_bounds = this->lower_ww_;
252✔
642
  auto& rel_err = this->upper_ww_;
252✔
643

644
  // get mesh volumes
645
  auto mesh_vols = this->mesh()->volumes();
252✔
646

647
  // Calculate mean (new_bounds) and relative error
648
#pragma omp parallel for collapse(2) schedule(static)
140✔
649
  for (int e = 0; e < e_bins; e++) {
934✔
650
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
651
      // Calculate mean
652
      new_bounds(e, m) = sum(e, m) / n;
1,189,249✔
653
      // Calculate relative error
654
      if (sum(e, m) > 0.0) {
1,189,249✔
655
        double mean_val = new_bounds(e, m);
101,480✔
656
        double variance = (sum_sq(e, m) / n - mean_val * mean_val) / (n - 1);
101,480✔
657
        rel_err(e, m) = std::sqrt(variance) / mean_val;
101,480✔
658
      } else {
659
        rel_err(e, m) = INFTY;
1,087,769✔
660
      }
661
      if (value == "rel_err") {
1,189,249✔
662
        new_bounds(e, m) = 1.0 / rel_err(e, m);
345,000✔
663
      }
664
    }
665
  }
666

667
  // Divide by volume of mesh elements
668
#pragma omp parallel for collapse(2) schedule(static)
140✔
669
  for (int e = 0; e < e_bins; e++) {
934✔
670
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
671
      new_bounds(e, m) /= mesh_vols[m];
1,189,249✔
672
    }
673
  }
674

675
  if (method == WeightWindowUpdateMethod::MAGIC) {
252✔
676
    // For MAGIC, weight windows are proportional to the forward fluxes.
677
    // We normalize weight windows independently for each energy group.
678

679
    // Find group maximum and normalize (per energy group)
680
    for (int e = 0; e < e_bins; e++) {
1,870✔
681
      double group_max = 0.0;
1,716✔
682

683
      // Find maximum value across all elements in this energy group
684
#pragma omp parallel for schedule(static) reduction(max : group_max)
936✔
685
      for (int64_t m = 0; m < mesh_bins; m++) {
1,092,505✔
686
        if (new_bounds(e, m) > group_max) {
1,091,725✔
687
          group_max = new_bounds(e, m);
2,520✔
688
        }
689
      }
690

691
      // Normalize values in this energy group by the maximum value
692
      if (group_max > 0.0) {
1,716✔
693
        double norm_factor = 1.0 / (2.0 * group_max);
1,683✔
694
#pragma omp parallel for schedule(static)
918✔
695
        for (int64_t m = 0; m < mesh_bins; m++) {
1,091,590✔
696
          new_bounds(e, m) *= norm_factor;
1,090,825✔
697
        }
698
      }
699
    }
700
  } else {
701
    // For FW-CADIS, weight windows are inversely proportional to the adjoint
702
    // fluxes. We normalize the weight windows across all energy groups.
703
#pragma omp parallel for collapse(2) schedule(static)
56✔
704
    for (int e = 0; e < e_bins; e++) {
84✔
705
      for (int64_t m = 0; m < mesh_bins; m++) {
97,566✔
706
        // Take the inverse, but are careful not to divide by zero
707
        if (new_bounds(e, m) != 0.0) {
97,524✔
708
          new_bounds(e, m) = 1.0 / new_bounds(e, m);
69,660✔
709
        } else {
710
          new_bounds(e, m) = 0.0;
27,864!
711
        }
712
      }
713
    }
714

715
    // Find the maximum value across all elements
716
    double max_val = 0.0;
98✔
717
#pragma omp parallel for collapse(2) schedule(static) reduction(max : max_val)
56✔
718
    for (int e = 0; e < e_bins; e++) {
84✔
719
      for (int64_t m = 0; m < mesh_bins; m++) {
97,566✔
720
        if (new_bounds(e, m) > max_val) {
97,524✔
721
          max_val = new_bounds(e, m);
405✔
722
        }
723
      }
724
    }
725

726
    // Parallel normalization
727
    if (max_val > 0.0) {
98✔
728
      double norm_factor = 1.0 / (2.0 * max_val);
68✔
729
#pragma omp parallel for collapse(2) schedule(static)
38✔
730
      for (int e = 0; e < e_bins; e++) {
60✔
731
        for (int64_t m = 0; m < mesh_bins; m++) {
69,690✔
732
          new_bounds(e, m) *= norm_factor;
69,660✔
733
        }
734
      }
735
    }
736
  }
737

738
  // Final processing
739
#pragma omp parallel for collapse(2) schedule(static)
140✔
740
  for (int e = 0; e < e_bins; e++) {
934✔
741
    for (int64_t m = 0; m < mesh_bins; m++) {
1,190,071✔
742
      // Values where the mean is zero should be ignored
743
      if (sum(e, m) <= 0.0) {
1,189,249✔
744
        new_bounds(e, m) = -1.0;
1,087,769✔
745
      }
746
      // Values where the relative error is higher than the threshold should be
747
      // ignored
748
      else if (rel_err(e, m) > threshold) {
101,480✔
749
        new_bounds(e, m) = -1.0;
1,420✔
750
      }
751
      // Set the upper bounds
752
      upper_ww_(e, m) = ratio * lower_ww_(e, m);
1,189,249✔
753
    }
754
  }
755
}
252✔
756

757
void WeightWindows::check_tally_update_compatibility(const Tally* tally)
252✔
758
{
759
  // define the set of allowed filters for the tally
760
  const std::set<FilterType> allowed_filters = {
761
    FilterType::MESH, FilterType::ENERGY, FilterType::PARTICLE};
252✔
762

763
  // retrieve a mapping of filter type to filter index for the tally
764
  auto filter_indices = tally->filter_indices();
252✔
765

766
  // a mesh filter is required for a tally used to update weight windows
767
  if (!filter_indices.count(FilterType::MESH)) {
252!
768
    fatal_error(
×
769
      "A mesh filter is required for a tally to update weight window bounds");
770
  }
771

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

775
  // make sure that all of the filters present on the tally are allowed
776
  for (auto filter_pair : filter_indices) {
964✔
777
    if (allowed_filters.find(filter_pair.first) == allowed_filters.end()) {
712!
778
      fatal_error(fmt::format("Invalid filter type '{}' found on tally "
×
779
                              "used for weight window generation.",
780
        model::tally_filters[tally->filters(filter_pair.second)]->type_str()));
×
781
    }
782
  }
783

784
  if (mesh_filter->mesh() != mesh_idx_) {
252!
785
    int32_t mesh_filter_id = model::meshes[mesh_filter->mesh()]->id();
×
786
    int32_t ww_mesh_id = model::meshes[this->mesh_idx_]->id();
×
787
    fatal_error(fmt::format("Mesh filter {} uses a different mesh ({}) than "
×
788
                            "weight window {} mesh ({})",
789
      mesh_filter->id(), mesh_filter_id, id_, ww_mesh_id));
×
790
  }
791

792
  // if an energy filter exists, make sure the energy grid matches that of this
793
  // weight window object
794
  if (auto energy_filter = tally->get_filter<EnergyFilter>()) {
252✔
795
    std::vector<double> filter_bins = energy_filter->bins();
230✔
796
    std::set<double> filter_e_bounds(
797
      energy_filter->bins().begin(), energy_filter->bins().end());
230✔
798
    if (filter_e_bounds.size() != energy_bounds().size()) {
230!
799
      fatal_error(
×
800
        fmt::format("Energy filter {} does not have the same number of energy "
×
801
                    "bounds ({}) as weight window object {} ({})",
802
          energy_filter->id(), filter_e_bounds.size(), id_,
×
803
          energy_bounds().size()));
×
804
    }
805

806
    for (auto e : energy_bounds()) {
2,252✔
807
      if (filter_e_bounds.count(e) == 0) {
2,022!
808
        fatal_error(fmt::format(
×
809
          "Energy bounds of filter {} and weight windows {} do not match",
810
          energy_filter->id(), id_));
×
811
      }
812
    }
813
  }
230✔
814
}
252✔
815

816
void WeightWindows::to_hdf5(hid_t group) const
134✔
817
{
818
  hid_t ww_group = create_group(group, fmt::format("weight_windows_{}", id()));
268✔
819

820
  write_dataset(ww_group, "mesh", this->mesh()->id());
134✔
821
  write_dataset(
134✔
822
    ww_group, "particle_type", openmc::particle_type_to_str(particle_type_));
268✔
823
  write_dataset(ww_group, "energy_bounds", energy_bounds_);
134✔
824
  write_dataset(ww_group, "lower_ww_bounds", lower_ww_);
134✔
825
  write_dataset(ww_group, "upper_ww_bounds", upper_ww_);
134✔
826
  write_dataset(ww_group, "survival_ratio", survival_ratio_);
134✔
827
  write_dataset(ww_group, "max_lower_bound_ratio", max_lb_ratio_);
134✔
828
  write_dataset(ww_group, "max_split", max_split_);
134✔
829
  write_dataset(ww_group, "weight_cutoff", weight_cutoff_);
134✔
830

831
  close_group(ww_group);
134✔
832
}
134✔
833

834
WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
82✔
835
{
836
  // read information from the XML node
837
  int32_t mesh_id = std::stoi(get_node_value(node, "mesh"));
82✔
838
  int32_t mesh_idx = model::mesh_map[mesh_id];
82✔
839
  max_realizations_ = std::stoi(get_node_value(node, "max_realizations"));
82✔
840

841
  int32_t active_batches = settings::n_batches - settings::n_inactive;
82✔
842
  if (max_realizations_ > active_batches) {
82✔
843
    auto msg =
844
      fmt::format("The maximum number of specified tally realizations ({}) is "
845
                  "greater than the number of active batches ({}).",
846
        max_realizations_, active_batches);
31✔
847
    warning(msg);
17✔
848
  }
17✔
849
  auto tmp_str = get_node_value(node, "particle_type", true, true);
82✔
850
  auto particle_type = str_to_particle_type(tmp_str);
82✔
851

852
  update_interval_ = std::stoi(get_node_value(node, "update_interval"));
82✔
853
  on_the_fly_ = get_node_value_bool(node, "on_the_fly");
82✔
854

855
  std::vector<double> e_bounds;
82✔
856
  if (check_for_node(node, "energy_bounds")) {
82✔
857
    e_bounds = get_node_array<double>(node, "energy_bounds");
23✔
858
  } else {
859
    int p_type = static_cast<int>(particle_type);
59✔
860
    e_bounds.push_back(data::energy_min[p_type]);
59✔
861
    e_bounds.push_back(data::energy_max[p_type]);
59✔
862
  }
863

864
  // set method
865
  std::string method_string = get_node_value(node, "method");
82✔
866
  if (method_string == "magic") {
82✔
867
    method_ = WeightWindowUpdateMethod::MAGIC;
33✔
868
    if (settings::solver_type == SolverType::RANDOM_RAY &&
33!
869
        FlatSourceDomain::adjoint_) {
870
      fatal_error("Random ray weight window generation with MAGIC cannot be "
×
871
                  "done in adjoint mode.");
872
    }
873
  } else if (method_string == "fw_cadis") {
49!
874
    method_ = WeightWindowUpdateMethod::FW_CADIS;
49✔
875
    if (settings::solver_type != SolverType::RANDOM_RAY) {
49!
876
      fatal_error("FW-CADIS can only be run in random ray solver mode.");
×
877
    }
878
    FlatSourceDomain::adjoint_ = true;
49✔
879
  } else {
880
    fatal_error(fmt::format(
×
881
      "Unknown weight window update method '{}' specified", method_string));
882
  }
883

884
  // parse non-default update parameters if specified
885
  if (check_for_node(node, "update_parameters")) {
82✔
886
    pugi::xml_node params_node = node.child("update_parameters");
22✔
887
    if (check_for_node(params_node, "value"))
22!
888
      tally_value_ = get_node_value(params_node, "value");
22✔
889
    if (check_for_node(params_node, "threshold"))
22!
890
      threshold_ = std::stod(get_node_value(params_node, "threshold"));
22✔
891
    if (check_for_node(params_node, "ratio")) {
22!
892
      ratio_ = std::stod(get_node_value(params_node, "ratio"));
22✔
893
    }
894
  }
895

896
  // check update parameter values
897
  if (tally_value_ != "mean" && tally_value_ != "rel_err") {
82!
898
    fatal_error(fmt::format("Unsupported tally value '{}' specified for "
×
899
                            "weight window generation.",
900
      tally_value_));
×
901
  }
902
  if (threshold_ <= 0.0)
82!
903
    fatal_error(fmt::format("Invalid relative error threshold '{}' (<= 0.0) "
×
904
                            "specified for weight window generation",
905
      ratio_));
×
906
  if (ratio_ <= 1.0)
82!
907
    fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) "
×
908
                            "specified for weight window generation"));
909

910
  // create a matching weight windows object
911
  auto wws = WeightWindows::create();
82✔
912
  ww_idx_ = wws->index();
82✔
913
  wws->set_mesh(mesh_idx);
82✔
914
  if (e_bounds.size() > 0)
82!
915
    wws->set_energy_bounds(e_bounds);
82✔
916
  wws->set_particle_type(particle_type);
82✔
917
  wws->set_defaults();
82✔
918
}
82✔
919

920
void WeightWindowsGenerator::create_tally()
82✔
921
{
922
  const auto& wws = variance_reduction::weight_windows[ww_idx_];
82✔
923

924
  // create a tally based on the WWG information
925
  Tally* ww_tally = Tally::create();
82✔
926
  tally_idx_ = model::tally_map[ww_tally->id()];
82✔
927
  ww_tally->set_scores({"flux"});
164!
928

929
  int32_t mesh_id = wws->mesh()->id();
82✔
930
  int32_t mesh_idx = model::mesh_map.at(mesh_id);
82✔
931
  // see if there's already a mesh filter using this mesh
932
  bool found_mesh_filter = false;
82✔
933
  for (const auto& f : model::tally_filters) {
259✔
934
    if (f->type() == FilterType::MESH) {
188✔
935
      const auto* mesh_filter = dynamic_cast<MeshFilter*>(f.get());
11!
936
      if (mesh_filter->mesh() == mesh_idx && !mesh_filter->translated() &&
22!
937
          !mesh_filter->rotated()) {
11!
938
        ww_tally->add_filter(f.get());
11✔
939
        found_mesh_filter = true;
11✔
940
        break;
11✔
941
      }
942
    }
943
  }
944

945
  if (!found_mesh_filter) {
82✔
946
    auto mesh_filter = Filter::create("mesh");
71✔
947
    openmc_mesh_filter_set_mesh(mesh_filter->index(), model::mesh_map[mesh_id]);
71✔
948
    ww_tally->add_filter(mesh_filter);
71✔
949
  }
950

951
  const auto& e_bounds = wws->energy_bounds();
82✔
952
  if (e_bounds.size() > 0) {
82!
953
    auto energy_filter = Filter::create("energy");
82✔
954
    openmc_energy_filter_set_bins(
164✔
955
      energy_filter->index(), e_bounds.size(), e_bounds.data());
82✔
956
    ww_tally->add_filter(energy_filter);
82✔
957
  }
958

959
  // add a particle filter
960
  auto particle_type = wws->particle_type();
82✔
961
  auto particle_filter = Filter::create("particle");
82✔
962
  auto pf = dynamic_cast<ParticleFilter*>(particle_filter);
82!
963
  pf->set_particles({&particle_type, 1});
82✔
964
  ww_tally->add_filter(particle_filter);
82✔
965
}
82✔
966

967
void WeightWindowsGenerator::update() const
2,417✔
968
{
969
  const auto& wws = variance_reduction::weight_windows[ww_idx_];
2,417✔
970

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

973
  // If in random ray mode, only update on the last batch
974
  if (settings::solver_type == SolverType::RANDOM_RAY) {
2,417✔
975
    if (simulation::current_batch != settings::n_batches) {
2,252✔
976
      return;
2,154✔
977
    }
978
    // If in Monte Carlo mode and beyond the number of max realizations or
979
    // not at the correct update interval, skip the update
980
  } else if (max_realizations_ < tally->n_realizations_ ||
165✔
981
             tally->n_realizations_ % update_interval_ != 0) {
33!
982
    return;
132✔
983
  }
984

985
  wws->update_weights(tally, tally_value_, threshold_, ratio_, method_);
131✔
986

987
  // if we're not doing on the fly generation, reset the tally results once
988
  // we're done with the update
989
  if (!on_the_fly_)
131!
UNCOV
990
    tally->reset();
×
991

992
  // TODO: deactivate or remove tally once weight window generation is
993
  // complete
994
}
995

996
//==============================================================================
997
// Non-member functions
998
//==============================================================================
999

1000
void finalize_variance_reduction()
7,624✔
1001
{
1002
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
7,706✔
1003
    wwg->create_tally();
82✔
1004
  }
1005
}
7,624✔
1006

1007
//==============================================================================
1008
// C API
1009
//==============================================================================
1010

1011
int verify_ww_index(int32_t index)
1,991✔
1012
{
1013
  if (index < 0 || index >= variance_reduction::weight_windows.size()) {
1,991!
1014
    set_errmsg(fmt::format("Index '{}' for weight windows is invalid", index));
×
UNCOV
1015
    return OPENMC_E_OUT_OF_BOUNDS;
×
1016
  }
1017
  return 0;
1,991✔
1018
}
1019

1020
extern "C" int openmc_get_weight_windows_index(int32_t id, int32_t* idx)
165✔
1021
{
1022
  auto it = variance_reduction::ww_map.find(id);
165✔
1023
  if (it == variance_reduction::ww_map.end()) {
165!
1024
    set_errmsg(fmt::format("No weight windows exist with ID={}", id));
×
UNCOV
1025
    return OPENMC_E_INVALID_ID;
×
1026
  }
1027

1028
  *idx = it->second;
165✔
1029
  return 0;
165✔
1030
}
1031

1032
extern "C" int openmc_weight_windows_get_id(int32_t index, int32_t* id)
517✔
1033
{
1034
  if (int err = verify_ww_index(index))
517!
UNCOV
1035
    return err;
×
1036

1037
  const auto& wws = variance_reduction::weight_windows.at(index);
517✔
1038
  *id = wws->id();
517✔
1039
  return 0;
517✔
1040
}
1041

1042
extern "C" int openmc_weight_windows_set_id(int32_t index, int32_t id)
154✔
1043
{
1044
  if (int err = verify_ww_index(index))
154!
UNCOV
1045
    return err;
×
1046

1047
  const auto& wws = variance_reduction::weight_windows.at(index);
154✔
1048
  wws->set_id(id);
154✔
1049
  return 0;
154✔
1050
}
1051

1052
extern "C" int openmc_weight_windows_update_magic(int32_t ww_idx,
121✔
1053
  int32_t tally_idx, const char* value, double threshold, double ratio)
1054
{
1055
  if (int err = verify_ww_index(ww_idx))
121!
UNCOV
1056
    return err;
×
1057

1058
  if (tally_idx < 0 || tally_idx >= model::tallies.size()) {
121!
1059
    set_errmsg(fmt::format("Index '{}' for tally is invalid", tally_idx));
×
UNCOV
1060
    return OPENMC_E_OUT_OF_BOUNDS;
×
1061
  }
1062

1063
  // get the requested tally
1064
  const Tally* tally = model::tallies.at(tally_idx).get();
121✔
1065

1066
  // get the WeightWindows object
1067
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
121✔
1068

1069
  wws->update_weights(tally, value, threshold, ratio);
121✔
1070

1071
  return 0;
121✔
1072
}
1073

1074
extern "C" int openmc_weight_windows_set_mesh(int32_t ww_idx, int32_t mesh_idx)
154✔
1075
{
1076
  if (int err = verify_ww_index(ww_idx))
154!
UNCOV
1077
    return err;
×
1078
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
154✔
1079
  wws->set_mesh(mesh_idx);
154✔
1080
  return 0;
154✔
1081
}
1082

1083
extern "C" int openmc_weight_windows_get_mesh(int32_t ww_idx, int32_t* mesh_idx)
11✔
1084
{
1085
  if (int err = verify_ww_index(ww_idx))
11!
UNCOV
1086
    return err;
×
1087
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
11✔
1088
  *mesh_idx = model::mesh_map.at(wws->mesh()->id());
11✔
1089
  return 0;
11✔
1090
}
1091

1092
extern "C" int openmc_weight_windows_set_energy_bounds(
132✔
1093
  int32_t ww_idx, double* e_bounds, size_t e_bounds_size)
1094
{
1095
  if (int err = verify_ww_index(ww_idx))
132!
UNCOV
1096
    return err;
×
1097
  const auto& wws = variance_reduction::weight_windows.at(ww_idx);
132✔
1098
  wws->set_energy_bounds({e_bounds, e_bounds_size});
132✔
1099
  return 0;
132✔
1100
}
1101

1102
extern "C" int openmc_weight_windows_get_energy_bounds(
11✔
1103
  int32_t ww_idx, const double** e_bounds, size_t* e_bounds_size)
1104
{
1105
  if (int err = verify_ww_index(ww_idx))
11!
UNCOV
1106
    return err;
×
1107
  const auto& wws = variance_reduction::weight_windows[ww_idx].get();
11✔
1108
  *e_bounds = wws->energy_bounds().data();
11✔
1109
  *e_bounds_size = wws->energy_bounds().size();
11✔
1110
  return 0;
11✔
1111
}
1112

1113
extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle)
176✔
1114
{
1115
  if (int err = verify_ww_index(index))
176!
UNCOV
1116
    return err;
×
1117

1118
  const auto& wws = variance_reduction::weight_windows.at(index);
176✔
1119
  wws->set_particle_type(static_cast<ParticleType>(particle));
176✔
1120
  return 0;
176✔
1121
}
1122

1123
extern "C" int openmc_weight_windows_get_particle(int32_t index, int* particle)
44✔
1124
{
1125
  if (int err = verify_ww_index(index))
44!
UNCOV
1126
    return err;
×
1127

1128
  const auto& wws = variance_reduction::weight_windows.at(index);
44✔
1129
  *particle = static_cast<int>(wws->particle_type());
44✔
1130
  return 0;
44✔
1131
}
1132

1133
extern "C" int openmc_weight_windows_get_bounds(int32_t index,
484✔
1134
  const double** lower_bounds, const double** upper_bounds, size_t* size)
1135
{
1136
  if (int err = verify_ww_index(index))
484!
UNCOV
1137
    return err;
×
1138

1139
  const auto& wws = variance_reduction::weight_windows[index];
484✔
1140
  *size = wws->lower_ww_bounds().size();
484✔
1141
  *lower_bounds = wws->lower_ww_bounds().data();
484✔
1142
  *upper_bounds = wws->upper_ww_bounds().data();
484✔
1143
  return 0;
484✔
1144
}
1145

1146
extern "C" int openmc_weight_windows_set_bounds(int32_t index,
11✔
1147
  const double* lower_bounds, const double* upper_bounds, size_t size)
1148
{
1149
  if (int err = verify_ww_index(index))
11!
UNCOV
1150
    return err;
×
1151

1152
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1153
  wws->set_bounds({lower_bounds, size}, {upper_bounds, size});
11✔
1154
  return 0;
11✔
1155
}
1156

1157
extern "C" int openmc_weight_windows_get_survival_ratio(
33✔
1158
  int32_t index, double* ratio)
1159
{
1160
  if (int err = verify_ww_index(index))
33!
UNCOV
1161
    return err;
×
1162
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1163
  *ratio = wws->survival_ratio();
33✔
1164
  return 0;
33✔
1165
}
1166

1167
extern "C" int openmc_weight_windows_set_survival_ratio(
11✔
1168
  int32_t index, double ratio)
1169
{
1170
  if (int err = verify_ww_index(index))
11!
UNCOV
1171
    return err;
×
1172
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1173
  wws->survival_ratio() = ratio;
11✔
1174
  std::cout << "Survival ratio: " << wws->survival_ratio() << std::endl;
11✔
1175
  return 0;
11✔
1176
}
1177

1178
extern "C" int openmc_weight_windows_get_max_lower_bound_ratio(
33✔
1179
  int32_t index, double* lb_ratio)
1180
{
1181
  if (int err = verify_ww_index(index))
33!
UNCOV
1182
    return err;
×
1183
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1184
  *lb_ratio = wws->max_lower_bound_ratio();
33✔
1185
  return 0;
33✔
1186
}
1187

1188
extern "C" int openmc_weight_windows_set_max_lower_bound_ratio(
11✔
1189
  int32_t index, double lb_ratio)
1190
{
1191
  if (int err = verify_ww_index(index))
11!
UNCOV
1192
    return err;
×
1193
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1194
  wws->max_lower_bound_ratio() = lb_ratio;
11✔
1195
  return 0;
11✔
1196
}
1197

1198
extern "C" int openmc_weight_windows_get_weight_cutoff(
33✔
1199
  int32_t index, double* cutoff)
1200
{
1201
  if (int err = verify_ww_index(index))
33!
UNCOV
1202
    return err;
×
1203
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1204
  *cutoff = wws->weight_cutoff();
33✔
1205
  return 0;
33✔
1206
}
1207

1208
extern "C" int openmc_weight_windows_set_weight_cutoff(
11✔
1209
  int32_t index, double cutoff)
1210
{
1211
  if (int err = verify_ww_index(index))
11!
UNCOV
1212
    return err;
×
1213
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1214
  wws->weight_cutoff() = cutoff;
11✔
1215
  return 0;
11✔
1216
}
1217

1218
extern "C" int openmc_weight_windows_get_max_split(
33✔
1219
  int32_t index, int* max_split)
1220
{
1221
  if (int err = verify_ww_index(index))
33!
UNCOV
1222
    return err;
×
1223
  const auto& wws = variance_reduction::weight_windows[index];
33✔
1224
  *max_split = wws->max_split();
33✔
1225
  return 0;
33✔
1226
}
1227

1228
extern "C" int openmc_weight_windows_set_max_split(int32_t index, int max_split)
11✔
1229
{
1230
  if (int err = verify_ww_index(index))
11!
UNCOV
1231
    return err;
×
1232
  const auto& wws = variance_reduction::weight_windows[index];
11✔
1233
  wws->max_split() = max_split;
11✔
1234
  return 0;
11✔
1235
}
1236

1237
extern "C" int openmc_extend_weight_windows(
154✔
1238
  int32_t n, int32_t* index_start, int32_t* index_end)
1239
{
1240
  if (index_start)
154!
1241
    *index_start = variance_reduction::weight_windows.size();
154✔
1242
  if (index_end)
154!
UNCOV
1243
    *index_end = variance_reduction::weight_windows.size() + n - 1;
×
1244
  for (int i = 0; i < n; ++i)
308✔
1245
    variance_reduction::weight_windows.push_back(make_unique<WeightWindows>());
154✔
1246
  return 0;
154✔
1247
}
1248

1249
extern "C" size_t openmc_weight_windows_size()
154✔
1250
{
1251
  return variance_reduction::weight_windows.size();
154✔
1252
}
1253

1254
extern "C" int openmc_weight_windows_export(const char* filename)
164✔
1255
{
1256

1257
  if (!mpi::master)
164✔
1258
    return 0;
30✔
1259

1260
  std::string name = filename ? filename : "weight_windows.h5";
268✔
1261

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

1264
  hid_t ww_file = file_open(name, 'w');
134✔
1265

1266
  // Write file type
1267
  write_attribute(ww_file, "filetype", "weight_windows");
134✔
1268

1269
  // Write revisiion number for state point file
1270
  write_attribute(ww_file, "version", VERSION_WEIGHT_WINDOWS);
134✔
1271

1272
  hid_t weight_windows_group = create_group(ww_file, "weight_windows");
134✔
1273

1274
  hid_t mesh_group = create_group(ww_file, "meshes");
134✔
1275

1276
  std::vector<int32_t> mesh_ids;
134✔
1277
  std::vector<int32_t> ww_ids;
134✔
1278
  for (const auto& ww : variance_reduction::weight_windows) {
268✔
1279

1280
    ww->to_hdf5(weight_windows_group);
134✔
1281
    ww_ids.push_back(ww->id());
134✔
1282

1283
    // if the mesh has already been written, move on
1284
    int32_t mesh_id = ww->mesh()->id();
134✔
1285
    if (std::find(mesh_ids.begin(), mesh_ids.end(), mesh_id) != mesh_ids.end())
134!
UNCOV
1286
      continue;
×
1287

1288
    mesh_ids.push_back(mesh_id);
134✔
1289
    ww->mesh()->to_hdf5(mesh_group);
134✔
1290
  }
1291

1292
  write_attribute(mesh_group, "n_meshes", mesh_ids.size());
134✔
1293
  write_attribute(mesh_group, "ids", mesh_ids);
134✔
1294
  close_group(mesh_group);
134✔
1295

1296
  write_attribute(weight_windows_group, "n_weight_windows", ww_ids.size());
134✔
1297
  write_attribute(weight_windows_group, "ids", ww_ids);
134✔
1298
  close_group(weight_windows_group);
134✔
1299

1300
  file_close(ww_file);
134✔
1301

1302
  return 0;
134✔
1303
}
134✔
1304

1305
extern "C" int openmc_weight_windows_import(const char* filename)
11✔
1306
{
1307
  std::string name = filename ? filename : "weight_windows.h5";
11!
1308

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

1312
  if (!file_exists(name)) {
11!
UNCOV
1313
    set_errmsg(fmt::format("File '{}' does not exist", name));
×
1314
  }
1315

1316
  hid_t ww_file = file_open(name, 'r');
11✔
1317

1318
  // Check that filetype is correct
1319
  std::string filetype;
11✔
1320
  read_attribute(ww_file, "filetype", filetype);
11✔
1321
  if (filetype != "weight_windows") {
11!
1322
    file_close(ww_file);
×
1323
    set_errmsg(fmt::format("File '{}' is not a weight windows file.", name));
×
UNCOV
1324
    return OPENMC_E_INVALID_ARGUMENT;
×
1325
  }
1326

1327
  // Check that the file version is compatible
1328
  std::array<int, 2> file_version;
1329
  read_attribute(ww_file, "version", file_version);
11✔
1330
  if (file_version[0] != VERSION_WEIGHT_WINDOWS[0]) {
11!
1331
    std::string err_msg =
1332
      fmt::format("File '{}' has version {} which is incompatible with the "
1333
                  "expected version ({}).",
1334
        name, file_version, VERSION_WEIGHT_WINDOWS);
×
1335
    set_errmsg(err_msg);
×
1336
    return OPENMC_E_INVALID_ARGUMENT;
×
UNCOV
1337
  }
×
1338

1339
  hid_t weight_windows_group = open_group(ww_file, "weight_windows");
11✔
1340

1341
  hid_t mesh_group = open_group(ww_file, "meshes");
11✔
1342

1343
  read_meshes(mesh_group);
11✔
1344

1345
  std::vector<std::string> names = group_names(weight_windows_group);
11✔
1346

1347
  for (const auto& name : names) {
22✔
1348
    WeightWindows::from_hdf5(weight_windows_group, name);
11✔
1349
  }
1350

1351
  close_group(weight_windows_group);
11✔
1352

1353
  file_close(ww_file);
11✔
1354

1355
  return 0;
11✔
1356
}
11✔
1357

1358
} // 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

© 2025 Coveralls, Inc