23084721708

Committed 14 Mar 2026 08:56AM UTC coverage: 81.599% (-0.5%) from 82.058%

Build # 23084721708

Build Type

Pull #2693

github

Committed by

web-flow

Commit Message

Merge 0ed23ee59 into bc9c31e0f

Pull Request Pull Request #2693: Add reactivity control to coupled transport-depletion analyses

Run Details

17575 of 25275 branches covered (69.54%)

Branch coverage included in aggregate %.

74 of 85 new or added lines in 4 files covered. (87.06%)

3755 existing lines in 99 files now uncovered.

58074 of 67433 relevant lines covered (86.12%)

47252067.37 hits per line

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

91.78

/src/volume_calc.cpp

#include "openmc/volume_calc.h"

#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/openmp_interface.h"
#include "openmc/output.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/timer.h"
#include "openmc/xml_interface.h"

#include "openmc/tensor.h"
#include <fmt/core.h>

#include <algorithm> // for copy
#include <cmath>     // for pow, sqrt
#include <unordered_set>

namespace openmc {

//==============================================================================
// Global variables
//==============================================================================

namespace model {
vector<VolumeCalculation> volume_calcs;
}

//==============================================================================
// VolumeCalculation implementation
//==============================================================================

VolumeCalculation::VolumeCalculation(pugi::xml_node node)
{
  // Read domain type (cell, material or universe)
  std::string domain_type = get_node_value(node, "domain_type");
  if (domain_type == "cell") {
    domain_type_ = TallyDomain::CELL;
  } else if (domain_type == "material") {
    domain_type_ = TallyDomain::MATERIAL;
  } else if (domain_type == "universe") {
    domain_type_ = TallyDomain::UNIVERSE;
  } else {
    fatal_error(std::string("Unrecognized domain type for stochastic "
                            "volume calculation: " +
                            domain_type));
  }

  // Read domain IDs, bounding corodinates and number of samples
  domain_ids_ = get_node_array<int>(node, "domain_ids");
  lower_left_ = get_node_array<double>(node, "lower_left");
  upper_right_ = get_node_array<double>(node, "upper_right");
  n_samples_ = std::stoull(get_node_value(node, "samples"));

  if (check_for_node(node, "threshold")) {
    pugi::xml_node threshold_node = node.child("threshold");

    threshold_ = std::stod(get_node_value(threshold_node, "threshold"));
    if (threshold_ <= 0.0) {
      fatal_error(fmt::format("Invalid error threshold {} provided for a "
                              "volume calculation.",
        threshold_));
    }

    std::string tmp = get_node_value(threshold_node, "type");
    if (tmp == "variance") {
      trigger_type_ = TriggerMetric::variance;
    } else if (tmp == "std_dev") {
      trigger_type_ = TriggerMetric::standard_deviation;
    } else if (tmp == "rel_err") {
      trigger_type_ = TriggerMetric::relative_error;
    } else {
      fatal_error(fmt::format(
        "Invalid volume calculation trigger type '{}' provided.", tmp));
    }
  }

  // Ensure there are no duplicates by copying elements to a set and then
  // comparing the length with the original vector
  std::unordered_set<int> unique_ids(domain_ids_.cbegin(), domain_ids_.cend());
  if (unique_ids.size() != domain_ids_.size()) {
    throw std::runtime_error {"Domain IDs for a volume calculation "
                              "must be unique."};
  }
}

vector<VolumeCalculation::Result> VolumeCalculation::execute() const
{
  // Check to make sure domain IDs are valid
  for (auto uid : domain_ids_) {
    switch (domain_type_) {
    case TallyDomain::CELL:
      if (model::cell_map.find(uid) == model::cell_map.end()) {
        throw std::runtime_error {fmt::format(
          "Cell {} in volume calculation does not exist in geometry.", uid)};
      }
      break;
    case TallyDomain::MATERIAL:
      if (model::material_map.find(uid) == model::material_map.end()) {
        throw std::runtime_error {fmt::format(
          "Material {} in volume calculation does not exist in geometry.",
          uid)};
      }
      break;
    case TallyDomain::UNIVERSE:
      if (model::universe_map.find(uid) == model::universe_map.end()) {
        throw std::runtime_error {fmt::format(
          "Universe {} in volume calculation does not exist in geometry.",
          uid)};
      }
    }
  }

  // Shared data that is collected from all threads
  int n = domain_ids_.size();
  vector<vector<uint64_t>> master_indices(
    n); // List of material indices for each domain
  vector<vector<uint64_t>> master_hits(
    n); // Number of hits for each material in each domain
  int iterations = 0;

  // Divide work over MPI processes
  uint64_t min_samples = n_samples_ / mpi::n_procs;
  uint64_t remainder = n_samples_ % mpi::n_procs;
  uint64_t i_start, i_end;
  if (mpi::rank < remainder) {
    i_start = (min_samples + 1) * mpi::rank;
    i_end = i_start + min_samples + 1;
  } else {
    i_start =
      (min_samples + 1) * remainder + (mpi::rank - remainder) * min_samples;
    i_end = i_start + min_samples;
  }

  while (true) {

#pragma omp parallel
    {
      // Variables that are private to each thread
      vector<vector<uint64_t>> indices(n);
      vector<vector<uint64_t>> hits(n);
      Particle p;

// Sample locations and count hits
#pragma omp for
      for (size_t i = i_start; i < i_end; i++) {
        uint64_t id = iterations * n_samples_ + i;
        uint64_t seed = init_seed(id, STREAM_VOLUME);

        p.n_coord() = 1;
        Position xi {prn(&seed), prn(&seed), prn(&seed)};
        p.r() = lower_left_ + xi * (upper_right_ - lower_left_);
        p.u() = {1. / std::sqrt(3.), 1. / std::sqrt(3.), 1. / std::sqrt(3.)};

        // If this location is not in the geometry at all, move on to next block
        if (!exhaustive_find_cell(p))
          continue;

        if (domain_type_ == TallyDomain::MATERIAL) {
          if (p.material() != MATERIAL_VOID) {
            for (int i_domain = 0; i_domain < n; i_domain++) {
              if (model::materials[p.material()]->id_ ==
                  domain_ids_[i_domain]) {
                this->check_hit(
                  p.material(), indices[i_domain], hits[i_domain]);
                break;
              }
            }
          }
        } else if (domain_type_ == TallyDomain::CELL) {
          for (int level = 0; level < p.n_coord(); ++level) {
            for (int i_domain = 0; i_domain < n; i_domain++) {
              if (model::cells[p.coord(level).cell()]->id_ ==
                  domain_ids_[i_domain]) {
                this->check_hit(
                  p.material(), indices[i_domain], hits[i_domain]);
                break;
              }
            }
          }
        } else if (domain_type_ == TallyDomain::UNIVERSE) {
          for (int level = 0; level < p.n_coord(); ++level) {
            for (int i_domain = 0; i_domain < n; ++i_domain) {
              if (model::universes[p.coord(level).universe()]->id_ ==
                  domain_ids_[i_domain]) {
                check_hit(p.material(), indices[i_domain], hits[i_domain]);
                break;
              }
            }
          }
        }
      }

      // At this point, each thread has its own pair of index/hits lists and we
      // now need to reduce them. OpenMP is not nearly smart enough to do this
      // on its own, so we have to manually reduce them
      for (int i_domain = 0; i_domain < n; ++i_domain) {
        reduce_indices_hits(indices[i_domain], hits[i_domain],
          master_indices[i_domain], master_hits[i_domain]);
      }
    } // omp parallel

    // Reduce hits onto master process

    // Determine volume of bounding box
    Position d {upper_right_ - lower_left_};
    double volume_sample = d.x * d.y * d.z;

    // bump iteration counter and get total number
    // of samples at this point
    iterations++;
    uint64_t total_samples = iterations * n_samples_;

    // warn user if total sample size is greater than what the uin64_t type can
    // represent
    if (total_samples == std::numeric_limits<uint64_t>::max()) {
      warning("The number of samples has exceeded the type used to track hits. "
              "Volume results may be inaccurate.");
    }

    // reset
    double trigger_val = -INFTY;

    // Set size for members of the Result struct
    vector<Result> results(n);

    for (int i_domain = 0; i_domain < n; ++i_domain) {
      // Get reference to result for this domain
      auto& result {results[i_domain]};

      // Create 2D array to store atoms/uncertainty for each nuclide. Later this
      // is compressed into vectors storing only those nuclides that are
      // non-zero
      auto n_nuc =
        settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();
      auto atoms =
        tensor::zeros<double>({static_cast<size_t>(n_nuc), size_t {2}});

#ifdef OPENMC_MPI
      if (mpi::master) {
        for (int j = 1; j < mpi::n_procs; j++) {
          int q;
          // retrieve results
          MPI_Recv(
            &q, 1, MPI_UINT64_T, j, 2 * j, mpi::intracomm, MPI_STATUS_IGNORE);
          vector<uint64_t> buffer(2 * q);
          MPI_Recv(buffer.data(), 2 * q, MPI_UINT64_T, j, 2 * j + 1,
            mpi::intracomm, MPI_STATUS_IGNORE);
          for (int k = 0; k < q; ++k) {
            bool already_added = false;
            for (int m = 0; m < master_indices[i_domain].size(); ++m) {
              if (buffer[2 * k] == master_indices[i_domain][m]) {
                master_hits[i_domain][m] += buffer[2 * k + 1];
                already_added = true;
                break;
              }
            }
            if (!already_added) {
              master_indices[i_domain].push_back(buffer[2 * k]);
              master_hits[i_domain].push_back(buffer[2 * k + 1]);
            }
          }
        }
      } else {
        int q = master_indices[i_domain].size();
        vector<uint64_t> buffer(2 * q);
        for (int k = 0; k < q; ++k) {
          buffer[2 * k] = master_indices[i_domain][k];
          buffer[2 * k + 1] = master_hits[i_domain][k];
        }

        MPI_Send(&q, 1, MPI_UINT64_T, 0, 2 * mpi::rank, mpi::intracomm);
        MPI_Send(buffer.data(), 2 * q, MPI_UINT64_T, 0, 2 * mpi::rank + 1,
          mpi::intracomm);
      }
#endif

      if (mpi::master) {
        size_t total_hits = 0;
        for (int j = 0; j < master_indices[i_domain].size(); ++j) {
          total_hits += master_hits[i_domain][j];
          double f =
            static_cast<double>(master_hits[i_domain][j]) / total_samples;
          double var_f = f * (1.0 - f) / total_samples;

          int i_material = master_indices[i_domain][j];
          if (i_material == MATERIAL_VOID)
            continue;

          const auto& mat = model::materials[i_material];
          for (int k = 0; k < mat->nuclide_.size(); ++k) {
            // Accumulate nuclide density
            int i_nuclide = mat->nuclide_[k];
            atoms(i_nuclide, 0) += mat->atom_density_[k] * f;
            atoms(i_nuclide, 1) += std::pow(mat->atom_density_[k], 2) * var_f;
          }
        }

        // Determine volume
        result.volume[0] =
          static_cast<double>(total_hits) / total_samples * volume_sample;
        result.volume[1] =
          std::sqrt(result.volume[0] * (volume_sample - result.volume[0]) /
                    total_samples);
        result.iterations = iterations;

        // update threshold value if needed
        if (trigger_type_ != TriggerMetric::not_active) {
          double val = 0.0;
          switch (trigger_type_) {
          case TriggerMetric::standard_deviation:
            val = result.volume[1];
            break;
          case TriggerMetric::relative_error:
            val = result.volume[0] == 0.0 ? INFTY
                                          : result.volume[1] / result.volume[0];
            break;
          case TriggerMetric::variance:
            val = result.volume[1] * result.volume[1];
            break;
          default:
            break;
          }
          // update max if entry is valid
          if (val > 0.0) {
            trigger_val = std::max(trigger_val, val);
          }
        }

        for (int j = 0; j < n_nuc; ++j) {
          // Determine total number of atoms. At this point, we have values in
          // atoms/b-cm. To get to atoms we multiply by 10^24 V.
          double mean = 1.0e24 * volume_sample * atoms(j, 0);
          double stdev = 1.0e24 * volume_sample * std::sqrt(atoms(j, 1));

          // Convert full arrays to vectors
          if (mean > 0.0) {
            result.nuclides.push_back(j);
            result.atoms.push_back(mean);
            result.uncertainty.push_back(stdev);
          }
        }
      }
    } // end domain loop

    // if no trigger is applied, we're done
    if (trigger_type_ == TriggerMetric::not_active) {
      return results;
    }

#ifdef OPENMC_MPI
    // update maximum error value on all processes
    MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm);
#endif

    // return results of the calculation
    if (trigger_val < threshold_) {
      return results;
    }

#ifdef OPENMC_MPI
    // if iterating in an MPI run, need to zero indices and hits so they aren't
    // counted twice
    if (!mpi::master) {
      for (auto& v : master_indices) {
        std::fill(v.begin(), v.end(), 0);
      }
      for (auto& v : master_hits) {
        std::fill(v.begin(), v.end(), 0);
      }
    }
#endif

  } // end while
}

void VolumeCalculation::to_hdf5(
  const std::string& filename, const vector<Result>& results) const
{
  // Create HDF5 file
  hid_t file_id = file_open(filename, 'w');

  // Write header info
  write_attribute(file_id, "filetype", "volume");
  write_attribute(file_id, "version", VERSION_VOLUME);
  write_attribute(file_id, "openmc_version", VERSION);
#ifdef GIT_SHA1
  write_attribute(file_id, "git_sha1", GIT_SHA1);
#endif

  // Write current date and time
  write_attribute(file_id, "date_and_time", time_stamp());

  // Write basic metadata
  write_attribute(file_id, "samples", n_samples_);
  write_attribute(file_id, "lower_left", lower_left_);
  write_attribute(file_id, "upper_right", upper_right_);
  // Write trigger info
  if (trigger_type_ != TriggerMetric::not_active) {
    write_attribute(file_id, "iterations", results[0].iterations);
    write_attribute(file_id, "threshold", threshold_);
    std::string trigger_str;
    switch (trigger_type_) {
    case TriggerMetric::variance:
      trigger_str = "variance";
      break;
    case TriggerMetric::standard_deviation:
      trigger_str = "std_dev";
      break;
    case TriggerMetric::relative_error:
      trigger_str = "rel_err";
      break;
    default:
      break;
    }
    write_attribute(file_id, "trigger_type", trigger_str);
  } else {
    write_attribute(file_id, "iterations", 1);
  }

  if (domain_type_ == TallyDomain::CELL) {
    write_attribute(file_id, "domain_type", "cell");
  } else if (domain_type_ == TallyDomain::MATERIAL) {
    write_attribute(file_id, "domain_type", "material");
  } else if (domain_type_ == TallyDomain::UNIVERSE) {
    write_attribute(file_id, "domain_type", "universe");
  }

  for (int i = 0; i < domain_ids_.size(); ++i) {
    hid_t group_id =
      create_group(file_id, fmt::format("domain_{}", domain_ids_[i]));

    // Write volume for domain
    const auto& result {results[i]};
    write_dataset(group_id, "volume", result.volume);

    // Create array of nuclide names from the vector
    auto n_nuc = result.nuclides.size();

    vector<std::string> nucnames;
    for (int i_nuc : result.nuclides) {
      nucnames.push_back(settings::run_CE ? data::nuclides[i_nuc]->name_
                                          : data::mg.nuclides_[i_nuc].name);
    }

    // Create array of total # of atoms with uncertainty for each nuclide
    tensor::Tensor<double> atom_data({static_cast<size_t>(n_nuc), size_t {2}});
    for (size_t k = 0; k < static_cast<size_t>(n_nuc); ++k) {
      atom_data(k, 0) = result.atoms[k];
      atom_data(k, 1) = result.uncertainty[k];
    }

    // Write results
    write_dataset(group_id, "nuclides", nucnames);
    write_dataset(group_id, "atoms", atom_data);

    close_group(group_id);
  }

  file_close(file_id);
}

void VolumeCalculation::check_hit(
  int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const
{

  // Check if this material was previously hit and if so, increment count
  bool already_hit = false;
  for (int j = 0; j < indices.size(); j++) {
    if (indices[j] == i_material) {
      hits[j]++;
      already_hit = true;
    }
  }

  // If the material was not previously hit, append an entry to the material
  // indices and hits lists
  if (!already_hit) {
    indices.push_back(i_material);
    hits.push_back(1);
  }
}

void free_memory_volume()
{
  openmc::model::volume_calcs.clear();
}

} // namespace openmc

//==============================================================================
// OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations
// that the user has specified and writes results to HDF5 files
//==============================================================================

int openmc_calculate_volumes()
{
  using namespace openmc;

  if (mpi::master) {
    header("STOCHASTIC VOLUME CALCULATION", 3);
  }
  Timer time_volume;
  time_volume.start();

  for (int i = 0; i < model::volume_calcs.size(); ++i) {
    write_message(4, "Running volume calculation {}", i + 1);

    // Run volume calculation
    const auto& vol_calc {model::volume_calcs[i]};
    std::vector<VolumeCalculation::Result> results;
    try {
      results = vol_calc.execute();
    } catch (const std::exception& e) {
      set_errmsg(e.what());
      return OPENMC_E_UNASSIGNED;
    }

    if (mpi::master) {
      std::string domain_type;
      if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {
        domain_type = "  Cell ";
      } else if (vol_calc.domain_type_ ==
                 VolumeCalculation::TallyDomain::MATERIAL) {
        domain_type = "  Material ";
      } else {
        domain_type = "  Universe ";
      }

      // Display domain volumes
      for (int j = 0; j < vol_calc.domain_ids_.size(); j++) {
        std::string region_name {""};
        if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {
          int cell_idx = model::cell_map[vol_calc.domain_ids_[j]];
          region_name = model::cells[cell_idx]->name();
        } else if (vol_calc.domain_type_ ==
                   VolumeCalculation::TallyDomain::MATERIAL) {
          int mat_idx = model::material_map[vol_calc.domain_ids_[j]];
          region_name = model::materials[mat_idx]->name();
        }
        if (region_name.size())
          region_name.insert(0, " "); // prepend space for formatting

        write_message(4, "{}{}{}: {} +/- {} cm^3", domain_type,
          vol_calc.domain_ids_[j], region_name, results[j].volume[0],
          results[j].volume[1]);
      }

      // Write volumes to HDF5 file
      std::string filename =
        fmt::format("{}volume_{}.h5", settings::path_output, i + 1);
      vol_calc.to_hdf5(filename, results);
    }
  }

  // Show elapsed time
  time_volume.stop();
  write_message(6, "Elapsed time: {} s", time_volume.elapsed());

  return 0;
}

1	#include "openmc/volume_calc.h"
2
3	#include "openmc/capi.h"
4	#include "openmc/cell.h"
5	#include "openmc/constants.h"
6	#include "openmc/error.h"
7	#include "openmc/geometry.h"
8	#include "openmc/hdf5_interface.h"
9	#include "openmc/material.h"
10	#include "openmc/message_passing.h"
11	#include "openmc/mgxs_interface.h"
12	#include "openmc/nuclide.h"
13	#include "openmc/openmp_interface.h"
14	#include "openmc/output.h"
15	#include "openmc/random_lcg.h"
16	#include "openmc/settings.h"
17	#include "openmc/timer.h"
18	#include "openmc/xml_interface.h"
19
20	#include "openmc/tensor.h"
21	#include <fmt/core.h>
22
23	#include <algorithm> // for copy
24	#include <cmath> // for pow, sqrt
25	#include <unordered_set>
26
27	namespace openmc {
28
29	//==============================================================================
30	// Global variables
31	//==============================================================================
32
33	namespace model {
34	vector<VolumeCalculation> volume_calcs;
35	}
36
37	//==============================================================================
38	// VolumeCalculation implementation
39	//==============================================================================
40
41	VolumeCalculation::VolumeCalculation(pugi::xml_node node)	314✔
42	{
43	// Read domain type (cell, material or universe)
44	std::string domain_type = get_node_value(node, "domain_type");	314✔
45	if (domain_type == "cell") {	314✔
46	domain_type_ = TallyDomain::CELL;	152✔
47	} else if (domain_type == "material") {	162✔
48	domain_type_ = TallyDomain::MATERIAL;	112✔
49	} else if (domain_type == "universe") {	50!
50	domain_type_ = TallyDomain::UNIVERSE;	50✔
51	} else {
UNCOV 52	fatal_error(std::string("Unrecognized domain type for stochastic "	×
53	"volume calculation: " +
54	domain_type));
55	}
56
57	// Read domain IDs, bounding corodinates and number of samples
58	domain_ids_ = get_node_array<int>(node, "domain_ids");	314✔
59	lower_left_ = get_node_array<double>(node, "lower_left");	314✔
60	upper_right_ = get_node_array<double>(node, "upper_right");	314✔
61	n_samples_ = std::stoull(get_node_value(node, "samples"));	628✔
62
63	if (check_for_node(node, "threshold")) {	314✔
64	pugi::xml_node threshold_node = node.child("threshold");	90✔
65
66	threshold_ = std::stod(get_node_value(threshold_node, "threshold"));	180✔
67	if (threshold_ <= 0.0) {	90!
UNCOV 68	fatal_error(fmt::format("Invalid error threshold {} provided for a "	×
69	"volume calculation.",
UNCOV 70	threshold_));	×
71	}
72
73	std::string tmp = get_node_value(threshold_node, "type");	90✔
74	if (tmp == "variance") {	90✔
75	trigger_type_ = TriggerMetric::variance;	30✔
76	} else if (tmp == "std_dev") {	60✔
77	trigger_type_ = TriggerMetric::standard_deviation;	30✔
78	} else if (tmp == "rel_err") {	30!
79	trigger_type_ = TriggerMetric::relative_error;	30✔
80	} else {
UNCOV 81	fatal_error(fmt::format(	×
82	"Invalid volume calculation trigger type '{}' provided.", tmp));
83	}
84	}	90✔
85
86	// Ensure there are no duplicates by copying elements to a set and then
87	// comparing the length with the original vector
88	std::unordered_set<int> unique_ids(domain_ids_.cbegin(), domain_ids_.cend());	314!
89	if (unique_ids.size() != domain_ids_.size()) {	156!
UNCOV 90	throw std::runtime_error {"Domain IDs for a volume calculation "	×
91	"must be unique."};	×
92	}
93	}	314✔
94
95	vector<VolumeCalculation::Result> VolumeCalculation::execute() const	300✔
96	{
97	// Check to make sure domain IDs are valid
98	for (auto uid : domain_ids_) {	846✔
99	switch (domain_type_) {	573!
100	case TallyDomain::CELL:	347✔
101	if (model::cell_map.find(uid) == model::cell_map.end()) {	347✔
102	throw std::runtime_error {fmt::format(	9!
103	"Cell {} in volume calculation does not exist in geometry.", uid)};	18✔
104	}
105	break;
106	case TallyDomain::MATERIAL:	176✔
107	if (model::material_map.find(uid) == model::material_map.end()) {	176✔
108	throw std::runtime_error {fmt::format(	9!
109	"Material {} in volume calculation does not exist in geometry.",
110	uid)};	18✔
111	}
112	break;
113	case TallyDomain::UNIVERSE:	50✔
114	if (model::universe_map.find(uid) == model::universe_map.end()) {	50✔
115	throw std::runtime_error {fmt::format(	9!
116	"Universe {} in volume calculation does not exist in geometry.",
117	uid)};	18✔
118	}
119	}
120	}
121
122	// Shared data that is collected from all threads
123	int n = domain_ids_.size();	273✔
124	vector<vector<uint64_t>> master_indices(	273✔
125	n); // List of material indices for each domain	546✔
126	vector<vector<uint64_t>> master_hits(	273✔
127	n); // Number of hits for each material in each domain	546✔
128	int iterations = 0;	273✔
129
130	// Divide work over MPI processes
131	uint64_t min_samples = n_samples_ / mpi::n_procs;	273✔
132	uint64_t remainder = n_samples_ % mpi::n_procs;	273✔
133	uint64_t i_start, i_end;	273✔
134	if (mpi::rank < remainder) {	273!
UNCOV 135	i_start = (min_samples + 1) * mpi::rank;	×
136	i_end = i_start + min_samples + 1;	×
137	} else {
138	i_start =	273✔
139	(min_samples + 1) * remainder + (mpi::rank - remainder) * min_samples;	273✔
140	i_end = i_start + min_samples;	273✔
141	}
142
143	while (true) {	19,293✔
144
145	#pragma omp parallel	11,573✔
146	{	7,720✔
147	// Variables that are private to each thread
148	vector<vector<uint64_t>> indices(n);	7,720✔
149	vector<vector<uint64_t>> hits(n);	7,720✔
150	Particle p;	7,720✔
151
152	// Sample locations and count hits
153	#pragma omp for
154	for (size_t i = i_start; i < i_end; i++) {	3,679,770✔
155	uint64_t id = iterations * n_samples_ + i;	3,672,050✔
156	uint64_t seed = init_seed(id, STREAM_VOLUME);	3,672,050✔
157
158	p.n_coord() = 1;	3,672,050✔
159	Position xi {prn(&seed), prn(&seed), prn(&seed)};	3,672,050✔
160	p.r() = lower_left_ + xi * (upper_right_ - lower_left_);	3,672,050✔
161	p.u() = {1. / std::sqrt(3.), 1. / std::sqrt(3.), 1. / std::sqrt(3.)};	3,672,050✔
162
163	// If this location is not in the geometry at all, move on to next block
164	if (!exhaustive_find_cell(p))	3,672,050✔
165	continue;	945,415✔
166
167	if (domain_type_ == TallyDomain::MATERIAL) {	2,726,635✔
168	if (p.material() != MATERIAL_VOID) {	761,515✔
169	for (int i_domain = 0; i_domain < n; i_domain++) {	1,422,380✔
170	if (model::materials[p.material()]->id_ ==	1,418,220✔
171	domain_ids_[i_domain]) {	1,418,220✔
172	this->check_hit(	755,090✔
173	p.material(), indices[i_domain], hits[i_domain]);	755,090✔
174	break;
175	}
176	}
177	}
178	} else if (domain_type_ == TallyDomain::CELL) {	1,965,120✔
179	for (int level = 0; level < p.n_coord(); ++level) {	2,443,400✔
180	for (int i_domain = 0; i_domain < n; i_domain++) {	1,456,285✔
181	if (model::cells[p.coord(level).cell()]->id_ ==	1,442,170✔
182	domain_ids_[i_domain]) {	1,442,170✔
183	this->check_hit(	1,209,090✔
184	p.material(), indices[i_domain], hits[i_domain]);	1,209,090✔
185	break;
186	}
187	}
188	}
189	} else if (domain_type_ == TallyDomain::UNIVERSE) {	744,925!
190	for (int level = 0; level < p.n_coord(); ++level) {	1,489,850✔
191	for (int i_domain = 0; i_domain < n; ++i_domain) {	744,925!
192	if (model::universes[p.coord(level).universe()]->id_ ==	744,925!
193	domain_ids_[i_domain]) {	744,925!
194	check_hit(p.material(), indices[i_domain], hits[i_domain]);	744,925✔
195	break;
196	}
197	}
198	}
199	}
200	}
201
202	// At this point, each thread has its own pair of index/hits lists and we
203	// now need to reduce them. OpenMP is not nearly smart enough to do this
204	// on its own, so we have to manually reduce them
205	for (int i_domain = 0; i_domain < n; ++i_domain) {	30,666✔
206	reduce_indices_hits(indices[i_domain], hits[i_domain],	22,946✔
207	master_indices[i_domain], master_hits[i_domain]);	22,946✔
208	}
209	} // omp parallel	7,720✔
210
211	// Reduce hits onto master process
212
213	// Determine volume of bounding box
214	Position d {upper_right_ - lower_left_};	19,293✔
215	double volume_sample = d.x * d.y * d.z;	19,293✔
216
217	// bump iteration counter and get total number
218	// of samples at this point
219	iterations++;	19,293✔
220	uint64_t total_samples = iterations * n_samples_;	19,293✔
221
222	// warn user if total sample size is greater than what the uin64_t type can
223	// represent
224	if (total_samples == std::numeric_limits<uint64_t>::max()) {	19,293!
UNCOV 225	warning("The number of samples has exceeded the type used to track hits. "	×
226	"Volume results may be inaccurate.");
227	}
228
229	// reset
230	double trigger_val = -INFTY;	19,293✔
231
232	// Set size for members of the Result struct
233	vector<Result> results(n);	19,293✔
234
235	for (int i_domain = 0; i_domain < n; ++i_domain) {	76,659✔
236	// Get reference to result for this domain
237	auto& result {results[i_domain]};	57,366✔
238
239	// Create 2D array to store atoms/uncertainty for each nuclide. Later this
240	// is compressed into vectors storing only those nuclides that are
241	// non-zero
242	auto n_nuc =	57,366✔
243	settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size();	57,366✔
244	auto atoms =	57,366✔
245	tensor::zeros<double>({static_cast<size_t>(n_nuc), size_t {2}});	57,366✔
246
247	#ifdef OPENMC_MPI
248	if (mpi::master) {	30,584✔
249	for (int j = 1; j < mpi::n_procs; j++) {	30,584✔
250	int q;	15,264✔
251	// retrieve results
252	MPI_Recv(	15,264✔
253	&q, 1, MPI_UINT64_T, j, 2 * j, mpi::intracomm, MPI_STATUS_IGNORE);
254	vector<uint64_t> buffer(2 * q);	15,264✔
255	MPI_Recv(buffer.data(), 2 * q, MPI_UINT64_T, j, 2 * j + 1,	15,264✔
256	mpi::intracomm, MPI_STATUS_IGNORE);
257	for (int k = 0; k < q; ++k) {	1,166,528✔
258	bool already_added = false;	2,287,264✔
259	for (int m = 0; m < master_indices[i_domain].size(); ++m) {	2,287,264✔
260	if (buffer[2 * k] == master_indices[i_domain][m]) {	2,287,224✔
261	master_hits[i_domain][m] += buffer[2 * k + 1];	1,151,224✔
262	already_added = true;	1,151,224✔
263	break;	1,151,224✔
264	}
265	}
266	if (!already_added) {	1,151,224✔
267	master_indices[i_domain].push_back(buffer[2 * k]);	40✔
268	master_hits[i_domain].push_back(buffer[2 * k + 1]);	40✔
269	}
270	}
271	}	15,264✔
272	} else {
273	int q = master_indices[i_domain].size();	15,264✔
274	vector<uint64_t> buffer(2 * q);	15,264✔
275	for (int k = 0; k < q; ++k) {	1,166,528✔
276	buffer[2 * k] = master_indices[i_domain][k];	1,151,264✔
277	buffer[2 * k + 1] = master_hits[i_domain][k];	1,151,264✔
278	}
279
280	MPI_Send(&q, 1, MPI_UINT64_T, 0, 2 * mpi::rank, mpi::intracomm);	15,264✔
281	MPI_Send(buffer.data(), 2 * q, MPI_UINT64_T, 0, 2 * mpi::rank + 1,	15,264✔
282	mpi::intracomm);
283	}	15,264✔
284	#endif
285
286	if (mpi::master) {	57,366✔
287	size_t total_hits = 0;
288	for (int j = 0; j < master_indices[i_domain].size(); ++j) {	89,290✔
289	total_hits += master_hits[i_domain][j];	47,188✔
290	double f =	47,188✔
291	static_cast<double>(master_hits[i_domain][j]) / total_samples;	47,188✔
292	double var_f = f * (1.0 - f) / total_samples;	47,188✔
293
294	int i_material = master_indices[i_domain][j];	47,188✔
295	if (i_material == MATERIAL_VOID)	47,188✔
296	continue;	11✔
297
298	const auto& mat = model::materials[i_material];	47,177✔
299	for (int k = 0; k < mat->nuclide_.size(); ++k) {	174,738✔
300	// Accumulate nuclide density
301	int i_nuclide = mat->nuclide_[k];	127,561✔
302	atoms(i_nuclide, 0) += mat->atom_density_[k] * f;	127,561✔
303	atoms(i_nuclide, 1) += std::pow(mat->atom_density_[k], 2) * var_f;	127,561✔
304	}
305	}
306
307	// Determine volume
308	result.volume[0] =	42,102✔
309	static_cast<double>(total_hits) / total_samples * volume_sample;	42,102✔
310	result.volume[1] =	42,102✔
311	std::sqrt(result.volume[0] * (volume_sample - result.volume[0]) /	42,102✔
312	total_samples);
313	result.iterations = iterations;	42,102✔
314
315	// update threshold value if needed
316	if (trigger_type_ != TriggerMetric::not_active) {	42,102✔
317	double val = 0.0;	41,844✔
318	switch (trigger_type_) {	41,844!
319	case TriggerMetric::standard_deviation:	34,518✔
320	val = result.volume[1];	34,518✔
321	break;	34,518✔
322	case TriggerMetric::relative_error:	396✔
323	val = result.volume[0] == 0.0 ? INFTY	396!
324	: result.volume[1] / result.volume[0];	396✔
325	break;	396✔
326	case TriggerMetric::variance:	6,930✔
327	val = result.volume[1] * result.volume[1];	6,930✔
328	break;	6,930✔
329	default:
330	break;
331	}
332	// update max if entry is valid
333	if (val > 0.0) {	41,844!
334	trigger_val = std::max(trigger_val, val);	55,858✔
335	}
336	}
337
338	for (int j = 0; j < n_nuc; ++j) {	253,184✔
339	// Determine total number of atoms. At this point, we have values in
340	// atoms/b-cm. To get to atoms we multiply by 10^24 V.
341	double mean = 1.0e24 * volume_sample * atoms(j, 0);	211,082✔
342	double stdev = 1.0e24 * volume_sample * std::sqrt(atoms(j, 1));	211,082✔
343
344	// Convert full arrays to vectors
345	if (mean > 0.0) {	211,082✔
346	result.nuclides.push_back(j);	112,345✔
347	result.atoms.push_back(mean);	112,345✔
348	result.uncertainty.push_back(stdev);	112,345✔
349	}
350	}
351	}
352	} // end domain loop	57,366✔
353
354	// if no trigger is applied, we're done
355	if (trigger_type_ == TriggerMetric::not_active) {	19,293✔
356	return results;	183✔
357	}
358
359	#ifdef OPENMC_MPI
360	// update maximum error value on all processes
361	MPI_Bcast(&trigger_val, 1, MPI_DOUBLE, 0, mpi::intracomm);	10,192✔
362	#endif
363
364	// return results of the calculation
365	if (trigger_val < threshold_) {	19,110✔
366	return results;	90✔
367	}
368
369	#ifdef OPENMC_MPI
370	// if iterating in an MPI run, need to zero indices and hits so they aren't
371	// counted twice
372	if (!mpi::master) {	10,144✔
373	for (auto& v : master_indices) {	20,224✔
374	std::fill(v.begin(), v.end(), 0);	30,304✔
375	}
376	for (auto& v : master_hits) {	20,224✔
377	std::fill(v.begin(), v.end(), 0);	30,304✔
378	}
379	}
380	#endif
381
382	} // end while	19,293✔
383	}
384
385	void VolumeCalculation::to_hdf5(	225✔
386	const std::string& filename, const vector<Result>& results) const
387	{
388	// Create HDF5 file
389	hid_t file_id = file_open(filename, 'w');	225✔
390
391	// Write header info
392	write_attribute(file_id, "filetype", "volume");	225✔
393	write_attribute(file_id, "version", VERSION_VOLUME);	225✔
394	write_attribute(file_id, "openmc_version", VERSION);	225✔
395	#ifdef GIT_SHA1
396	write_attribute(file_id, "git_sha1", GIT_SHA1);
397	#endif
398
399	// Write current date and time
400	write_attribute(file_id, "date_and_time", time_stamp());	450✔
401
402	// Write basic metadata
403	write_attribute(file_id, "samples", n_samples_);	225✔
404	write_attribute(file_id, "lower_left", lower_left_);	225✔
405	write_attribute(file_id, "upper_right", upper_right_);	225✔
406	// Write trigger info
407	if (trigger_type_ != TriggerMetric::not_active) {	225✔
408	write_attribute(file_id, "iterations", results[0].iterations);	66✔
409	write_attribute(file_id, "threshold", threshold_);	66✔
410	std::string trigger_str;	66!
411	switch (trigger_type_) {	66!
412	case TriggerMetric::variance:	22✔
413	trigger_str = "variance";	22✔
414	break;
415	case TriggerMetric::standard_deviation:	22✔
416	trigger_str = "std_dev";	22✔
417	break;
418	case TriggerMetric::relative_error:	22✔
419	trigger_str = "rel_err";	22✔
420	break;
421	default:
422	break;
423	}
424	write_attribute(file_id, "trigger_type", trigger_str);	132✔
425	} else {	66✔
426	write_attribute(file_id, "iterations", 1);	159✔
427	}
428
429	if (domain_type_ == TallyDomain::CELL) {	225✔
430	write_attribute(file_id, "domain_type", "cell");	112✔
431	} else if (domain_type_ == TallyDomain::MATERIAL) {	113✔
432	write_attribute(file_id, "domain_type", "material");	80✔
433	} else if (domain_type_ == TallyDomain::UNIVERSE) {	33!
434	write_attribute(file_id, "domain_type", "universe");	33✔
435	}
436
437	for (int i = 0; i < domain_ids_.size(); ++i) {	659✔
438	hid_t group_id =	434✔
439	create_group(file_id, fmt::format("domain_{}", domain_ids_[i]));	434✔
440
441	// Write volume for domain
442	const auto& result {results[i]};	434✔
443	write_dataset(group_id, "volume", result.volume);	434✔
444
445	// Create array of nuclide names from the vector
446	auto n_nuc = result.nuclides.size();	434✔
447
448	vector<std::string> nucnames;	434✔
449	for (int i_nuc : result.nuclides) {	1,723✔
450	nucnames.push_back(settings::run_CE ? data::nuclides[i_nuc]->name_	1,289✔
451	: data::mg.nuclides_[i_nuc].name);	429✔
452	}
453
454	// Create array of total # of atoms with uncertainty for each nuclide
455	tensor::Tensor<double> atom_data({static_cast<size_t>(n_nuc), size_t {2}});	434✔
456	for (size_t k = 0; k < static_cast<size_t>(n_nuc); ++k) {	1,723✔
457	atom_data(k, 0) = result.atoms[k];	1,289✔
458	atom_data(k, 1) = result.uncertainty[k];	1,289✔
459	}
460
461	// Write results
462	write_dataset(group_id, "nuclides", nucnames);	434✔
463	write_dataset(group_id, "atoms", atom_data);	434✔
464
465	close_group(group_id);	434✔
466	}	434✔
467
468	file_close(file_id);	225✔
469	}	225✔
470
471	void VolumeCalculation::check_hit(	10,960,031✔
472	int i_material, vector<uint64_t>& indices, vector<uint64_t>& hits) const
473	{
474
475	// Check if this material was previously hit and if so, increment count
476	bool already_hit = false;	10,960,031✔
477	for (int j = 0; j < indices.size(); j++) {	23,478,229✔
478	if (indices[j] == i_material) {	12,518,198✔
479	hits[j]++;	10,879,627✔
480	already_hit = true;	10,879,627✔
481	}
482	}
483
484	// If the material was not previously hit, append an entry to the material
485	// indices and hits lists
486	if (!already_hit) {	10,960,031✔
487	indices.push_back(i_material);	80,404✔
488	hits.push_back(1);	80,404✔
489	}
490	}	10,960,031✔
491
492	void free_memory_volume()	8,342✔
493	{
494	openmc::model::volume_calcs.clear();	8,342✔
495	}	8,342✔
496
497	} // namespace openmc
498
499	//==============================================================================
500	// OPENMC_CALCULATE_VOLUMES runs each of the stochastic volume calculations
501	// that the user has specified and writes results to HDF5 files
502	//==============================================================================
503
504	int openmc_calculate_volumes()	126✔
505	{
506	using namespace openmc;	126✔
507
508	if (mpi::master) {	126✔
509	header("STOCHASTIC VOLUME CALCULATION", 3);	118✔
510	}
511	Timer time_volume;	126✔
512	time_volume.start();	126✔
513
514	for (int i = 0; i < model::volume_calcs.size(); ++i) {	399✔
515	write_message(4, "Running volume calculation {}", i + 1);	300✔
516
517	// Run volume calculation
518	const auto& vol_calc {model::volume_calcs[i]};	300✔
519	std::vector<VolumeCalculation::Result> results;	300✔
520	try {	300✔
521	results = vol_calc.execute();	573✔
522	} catch (const std::exception& e) {	27!
523	set_errmsg(e.what());	27✔
524	return OPENMC_E_UNASSIGNED;	27✔
525	}	27✔
526
527	if (mpi::master) {	273✔
528	std::string domain_type;	225✔
529	if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {	225✔
530	domain_type = " Cell ";	112✔
531	} else if (vol_calc.domain_type_ ==	113✔
532	VolumeCalculation::TallyDomain::MATERIAL) {
533	domain_type = " Material ";	80✔
534	} else {
535	domain_type = " Universe ";	33✔
536	}
537
538	// Display domain volumes
539	for (int j = 0; j < vol_calc.domain_ids_.size(); j++) {	659✔
540	std::string region_name {""};	434✔
541	if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {	434✔
542	int cell_idx = model::cell_map[vol_calc.domain_ids_[j]];	266✔
543	region_name = model::cells[cell_idx]->name();	266✔
544	} else if (vol_calc.domain_type_ ==	168✔
545	VolumeCalculation::TallyDomain::MATERIAL) {
546	int mat_idx = model::material_map[vol_calc.domain_ids_[j]];	135✔
547	region_name = model::materials[mat_idx]->name();	135✔
548	}
549	if (region_name.size())	434✔
550	region_name.insert(0, " "); // prepend space for formatting	80✔
551
552	write_message(4, "{}{}{}: {} +/- {} cm^3", domain_type,	868✔
553	vol_calc.domain_ids_[j], region_name, results[j].volume[0],	434✔
554	results[j].volume[1]);	434✔
555	}	434✔
556
557	// Write volumes to HDF5 file
558	std::string filename =	225✔
559	fmt::format("{}volume_{}.h5", settings::path_output, i + 1);	225✔
560	vol_calc.to_hdf5(filename, results);	225✔
561	}	225✔
562	}	300✔
563
564	// Show elapsed time
565	time_volume.stop();	99✔
566	write_message(6, "Elapsed time: {} s", time_volume.elapsed());	99✔
567
568	return 0;	99✔
569	}

openmc-dev / openmc / 23084721708

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

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