refactor particle source import/export

This commit is contained in:
Gavin Ridley 2023-03-31 13:46:42 -04:00
parent a4f498c54f
commit 2172907bcf
13 changed files with 248 additions and 126 deletions

View file

@ -1,6 +1,8 @@
#ifndef OPENMC_FILE_UTILS_H
#define OPENMC_FILE_UTILS_H
#include <algorithm> // any_of
#include <cctype> // for isalpha
#include <fstream> // for ifstream
#include <string>
#include <sys/stat.h>
@ -33,6 +35,30 @@ inline bool file_exists(const std::string& filename)
return s.good();
}
// Gets the file extension of whatever string is passed in. This is defined as
// a sequence of strictly alphanumeric characters which follow the last period,
// i.e. at least one alphabet character is present, and zero or more numbers.
// If such a sequence of characters is not found, an empty string is returned.
inline std::string get_file_extension(const std::string& filename)
{
// check that at least one letter is present
const std::string::size_type last_period_pos = filename.find_last_of('.');
// no file extension
if (last_period_pos == std::string::npos)
return "";
const std::string ending = filename.substr(last_period_pos + 1);
// check that at least one character is present.
const bool has_alpha = std::any_of(ending.begin(), ending.end(),
[](char x) { return static_cast<bool>(std::isalpha(x)); });
if (has_alpha)
return ending;
else
return "";
}
} // namespace openmc
#endif // OPENMC_FILE_UTILS_H

View file

@ -4,6 +4,8 @@
#include "openmc/particle_data.h"
#include "openmc/vector.h"
#include <gsl/gsl-lite.hpp>
#include <string>
namespace openmc {
@ -26,11 +28,17 @@ vector<SourceSite> mcpl_source_sites(std::string path);
//! Write an MCPL source file
//
//! \param[in] filename Path to MCPL file
//! \param[in] surf_source_bank Whether to use the surface source bank
void write_mcpl_source_point(
const char* filename, bool surf_source_bank = false);
//! \param[in] filename Path to MCPL file
//! \param[in] source_bank Vector of SourceSites to write to file for this
//! MPI rank
//! \param[in] bank_indx Pointer to vector of site index ranges over all
//! MPI ranks, allowed to leave this argument alone
//! or pass a nullptr if not running in MPI mode.
//! The option of passing a nullptr has been left
//! for developers to experiment with serial code
//! before writing the fully MPI-parallel version.
void write_mcpl_source_point(const char* filename,
gsl::span<SourceSite> source_bank, vector<int64_t> const& bank_index);
} // namespace openmc
#endif // OPENMC_MCPL_INTERFACE_H

View file

@ -1,10 +1,14 @@
#ifndef OPENMC_MESSAGE_PASSING_H
#define OPENMC_MESSAGE_PASSING_H
#include <cstdint>
#ifdef OPENMC_MPI
#include <mpi.h>
#endif
#include "openmc/vector.h"
namespace openmc {
namespace mpi {
@ -17,6 +21,18 @@ extern MPI_Datatype source_site;
extern MPI_Comm intracomm;
#endif
// Calculates global indices of the bank particles
// across all ranks using a parallel scan. This is used to write
// the surface source file in parallel runs. It will probably
// be used in the future for other types of bank like particles
// in flight used to kick off transient simulations.
//
// More abstractly, this just takes a number from each MPI rank,
// and returns a vector which is the exclusive parallel scan across
// all of those numbers, having a length of the number of MPI ranks
// plus one.
vector<int64_t> calculate_parallel_index_vector(const int64_t size);
} // namespace mpi
} // namespace openmc

View file

@ -115,6 +115,12 @@ public:
T* data() { return data_.get(); }
const T* data() const { return data_.get(); }
//! Classic iterators
T* begin() { return data_.get(); }
const T* cbegin() const { return data_.get(); }
T* end() { return data_.get() + size_; }
const T* cend() const { return data_.get() + size_; }
private:
//==========================================================================
// Data members

View file

@ -3,18 +3,39 @@
#include <cstdint>
#include <gsl/gsl-lite.hpp>
#include "hdf5.h"
#include "openmc/capi.h"
#include "openmc/particle.h"
#include "openmc/shared_array.h"
#include "openmc/vector.h"
namespace openmc {
void load_state_point();
vector<int64_t> calculate_surf_source_size();
void write_source_point(const char* filename, bool surf_source_bank = false);
void write_source_bank(hid_t group_id, bool surf_source_bank);
// By passing in a filename, source bank, and list of source indices
// on each MPI rank, this writes an HDF5 file which contains that
// information which can later be read in by read_source_bank
// (defined below). If you're writing code to write out a new kind
// of particle bank, this function is the one you want to use!
//
// For example, this is used to write both the surface source sites
// or fission source sites for eigenvalue continuation runs.
//
// This function ends up calling write_source_bank, and is responsible
// for opening the file to be written to and controlling whether the
// write is done in parallel (if compiled with parallel HDF5).
void write_source_point(const char* filename, gsl::span<SourceSite> source_bank,
vector<int64_t> const& bank_index);
// This appends a source bank specification to an HDF5 file
// that's already open. It is used internally by write_source_point.
void write_source_bank(hid_t group_id, gsl::span<SourceSite> source_bank,
vector<int64_t> const& bank_index);
void read_source_bank(
hid_t group_id, vector<SourceSite>& sites, bool distribute);
void write_tally_results_nr(hid_t file_id);

View file

@ -2,6 +2,7 @@
#include "openmc/bank.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
@ -110,37 +111,18 @@ vector<SourceSite> mcpl_source_sites(std::string path)
//==============================================================================
#ifdef OPENMC_MCPL
void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
void write_mcpl_source_bank(mcpl_outfile_t file_id,
gsl::span<SourceSite> source_bank, vector<int64_t> const& bank_index)
{
int64_t dims_size = settings::n_particles;
int64_t count_size = simulation::work_per_rank;
// Set vectors for source bank and starting bank index of each process
vector<int64_t>* bank_index = &simulation::work_index;
vector<SourceSite>* source_bank = &simulation::source_bank;
vector<int64_t> surf_source_index_vector;
vector<SourceSite> surf_source_bank_vector;
if (surf_source_bank) {
surf_source_index_vector = calculate_surf_source_size();
dims_size = surf_source_index_vector[mpi::n_procs];
count_size = simulation::surf_source_bank.size();
bank_index = &surf_source_index_vector;
// Copy data in a SharedArray into a vector.
surf_source_bank_vector.resize(count_size);
surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
simulation::surf_source_bank.data() + count_size);
source_bank = &surf_source_bank_vector;
}
if (mpi::master) {
// Particles are writeen to disk from the master node only
// Save source bank sites since the array is overwritten below
#ifdef OPENMC_MPI
vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
vector<SourceSite> temp_source {source_bank.begin(), source_bank.end()};
#endif
// loop over the other nodes and receive data - then write those.
@ -151,11 +133,11 @@ void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
#ifdef OPENMC_MPI
if (i > 0)
MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
MPI_Recv(source_bank.data(), count[0], mpi::source_site, i, i,
mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// now write the source_bank data again.
for (const auto& site : *source_bank) {
for (const auto& site : source_bank) {
// particle is now at the iterator
// write it to the mcpl-file
mcpl_particle_t p;
@ -194,11 +176,11 @@ void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
}
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
std::copy(temp_source.begin(), temp_source.end(), source_bank.begin());
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
MPI_Send(source_bank.data(), count_size, mpi::source_site, 0, mpi::rank,
mpi::intracomm);
#endif
}
@ -207,17 +189,16 @@ void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
//==============================================================================
void write_mcpl_source_point(const char* filename, bool surf_source_bank)
void write_mcpl_source_point(const char* filename,
gsl::span<SourceSite> source_bank, vector<int64_t> const& bank_index)
{
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output,
simulation::current_batch, w);
std::string filename_(filename);
const auto extension = get_file_extension(filename_);
if (extension == "") {
filename_.append(".mcpl");
} else if (extension != "mcpl") {
fatal_error("write_mcpl_source_point was passed an incorrect file "
"extension. Must either pass nothing or .mcpl");
}
#ifdef OPENMC_MCPL
@ -236,7 +217,7 @@ void write_mcpl_source_point(const char* filename, bool surf_source_bank)
mcpl_hdr_set_srcname(file_id, line.c_str());
}
write_mcpl_source_bank(file_id, surf_source_bank);
write_mcpl_source_bank(file_id, source_bank, bank_index);
if (mpi::master) {
mcpl_close_outfile(file_id);

View file

@ -17,6 +17,29 @@ extern "C" bool openmc_master()
return mpi::master;
}
vector<int64_t> calculate_parallel_index_vector(const int64_t size)
{
vector<int64_t> result;
result.reserve(n_procs + 1);
#ifdef OPENMC_MPI
result.resize(n_procs);
vector<int64_t> bank_size(n_procs);
// Populate the result with cumulative sum of the number of
// surface source banks per process
MPI_Scan(&size, bank_size.data(), 1, MPI_INT64_T, MPI_SUM, intracomm);
MPI_Allgather(
bank_size.data(), 1, MPI_INT64_T, result.data(), 1, MPI_INT64_T, intracomm);
result.insert(result.begin(), 0);
#else
result.push_back(0);
result.push_back(size);
#endif
return result;
}
} // namespace mpi
} // namespace openmc

View file

@ -360,6 +360,8 @@ void print_build_info()
std::string png(n);
std::string profiling(n);
std::string coverage(n);
std::string mcpl(n);
std::string ncrystal(n);
#ifdef PHDF5
phdf5 = y;
@ -382,6 +384,12 @@ void print_build_info()
#ifdef COVERAGEBUILD
coverage = y;
#endif
#ifdef NCRYSTAL
ncrystal = y;
#endif
#ifdef OPENMC_MCPL
mcpl = y;
#endif
// Wraps macro variables in quotes
#define STRINGIFY(x) STRINGIFY2(x)
@ -396,6 +404,8 @@ void print_build_info()
fmt::print("PNG support: {}\n", png);
fmt::print("DAGMC support: {}\n", dagmc);
fmt::print("libMesh support: {}\n", libmesh);
fmt::print("NCrystal support: {}\n", profiling);
fmt::print("MCPL support: {}\n", profiling);
fmt::print("Coverage testing: {}\n", coverage);
fmt::print("Profiling flags: {}\n", profiling);
}

View file

@ -392,21 +392,32 @@ void finalize_batch()
// Write out a separate source point if it's been specified for this batch
if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
settings::source_write && settings::source_separate) {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
settings::path_output, simulation::current_batch, w);
gsl::span<SourceSite> bankspan(simulation::source_bank);
if (settings::source_mcpl_write) {
write_mcpl_source_point(nullptr);
write_mcpl_source_point(
source_point_filename.c_str(), bankspan, simulation::work_index);
} else {
write_source_point(nullptr);
write_source_point(
source_point_filename.c_str(), bankspan, simulation::work_index);
}
}
// Write a continously-overwritten source point if requested.
if (settings::source_latest) {
// note: correct file extension appended automatically
auto filename = settings::path_output + "source";
gsl::span<SourceSite> bankspan(simulation::source_bank);
if (settings::source_mcpl_write) {
auto filename = settings::path_output + "source.mcpl";
write_mcpl_source_point(filename.c_str());
write_mcpl_source_point(
filename.c_str(), bankspan, simulation::work_index);
} else {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
write_source_point(filename.c_str(), bankspan, simulation::work_index);
}
}
}
@ -414,12 +425,15 @@ void finalize_batch()
// Write out surface source if requested.
if (settings::surf_source_write &&
simulation::current_batch == settings::n_batches) {
auto filename = settings::path_output + "surface_source";
auto surf_work_index =
mpi::calculate_parallel_index_vector(simulation::surf_source_bank.size());
gsl::span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
simulation::surf_source_bank.size());
if (settings::surf_mcpl_write) {
auto filename = settings::path_output + "surface_source.mcpl";
write_mcpl_source_point(filename.c_str(), true);
write_mcpl_source_point(filename.c_str(), surfbankspan, surf_work_index);
} else {
auto filename = settings::path_output + "surface_source.h5";
write_source_point(filename.c_str(), true);
write_source_point(filename.c_str(), surfbankspan, surf_work_index);
}
}
}

View file

@ -403,7 +403,7 @@ void initialize_source()
write_message("Writing out initial source...", 5);
std::string filename = settings::path_output + "initial_source.h5";
hid_t file_id = file_open(filename, 'w', true);
write_source_bank(file_id, false);
write_source_bank(file_id, simulation::source_bank, simulation::work_index);
file_close(file_id);
}
}

View file

@ -13,6 +13,7 @@
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
@ -34,7 +35,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
{
simulation::time_statepoint.start();
// Set the filename
// If a nullptr is passed in, we assume that the user
// wants a default name for this, of the form like output/statepoint.20.h5
std::string filename_;
if (filename) {
filename_ = filename;
@ -47,6 +49,15 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
simulation::current_batch, w);
}
// If a file name was specified, ensure it has .h5 file extension
const auto extension = get_file_extension(filename_);
if (extension == "") {
filename_.append(".h5");
} else if (extension != "h5") {
fatal_error("openmc_statepoint_write was passed an incorrect file "
"extension. Must either have no file extension or .h5");
}
// Determine whether or not to write the source bank
bool write_source_ = write_source ? *write_source : true;
@ -316,7 +327,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
if (write_source_) {
if (mpi::master || parallel)
file_id = file_open(filename_, 'a', true);
write_source_bank(file_id, false);
write_source_bank(file_id, simulation::source_bank, simulation::work_index);
if (mpi::master || parallel)
file_close(file_id);
}
@ -541,31 +552,8 @@ hid_t h5banktype()
return banktype;
}
vector<int64_t> calculate_surf_source_size()
{
vector<int64_t> surf_source_index;
surf_source_index.reserve(mpi::n_procs + 1);
#ifdef OPENMC_MPI
surf_source_index.resize(mpi::n_procs);
vector<int64_t> bank_size(mpi::n_procs);
// Populate the surf_source_index with cumulative sum of the number of
// surface source banks per process
int64_t size = simulation::surf_source_bank.size();
MPI_Scan(&size, bank_size.data(), 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
MPI_Allgather(bank_size.data(), 1, MPI_INT64_T, surf_source_index.data(), 1,
MPI_INT64_T, mpi::intracomm);
surf_source_index.insert(surf_source_index.begin(), 0);
#else
surf_source_index.push_back(0);
surf_source_index.push_back(simulation::surf_source_bank.size());
#endif
return surf_source_index;
}
void write_source_point(const char* filename, bool surf_source_bank)
void write_source_point(const char* filename, gsl::span<SourceSite> source_bank,
vector<int64_t> const& bank_index)
{
// When using parallel HDF5, the file is written to collectively by all
// processes. With MPI-only, the file is opened and written by the master
@ -577,58 +565,59 @@ void write_source_point(const char* filename, bool surf_source_bank)
bool parallel = false;
#endif
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
if (!filename)
fatal_error("write_source_point filename needs a nonempty name.");
filename_ = fmt::format("{0}source.{1:0{2}}.h5", settings::path_output,
simulation::current_batch, w);
std::string filename_(filename);
const auto extension = get_file_extension(filename_);
if (extension == "") {
filename_.append(".h5");
} else if (extension != "h5") {
fatal_error("openmc_source_point was passed an incorrect file "
"extension. Must either have no file extension or .h5");
}
hid_t file_id;
if (mpi::master || parallel) {
file_id = file_open(filename_, 'w', true);
file_id = file_open(filename_.c_str(), 'w', true);
write_attribute(file_id, "filetype", "source");
}
// Get pointer to source bank and write to file
write_source_bank(file_id, surf_source_bank);
write_source_bank(file_id, source_bank, bank_index);
if (mpi::master || parallel)
file_close(file_id);
}
void write_source_bank(hid_t group_id, bool surf_source_bank)
void write_source_bank(hid_t group_id, gsl::span<SourceSite> source_bank,
vector<int64_t> const& bank_index)
{
hid_t banktype = h5banktype();
// Set total and individual process dataspace sizes for source bank
int64_t dims_size = settings::n_particles;
int64_t count_size = simulation::work_per_rank;
// TODO: are these correct?
// the old code was
//
// int64_t dims_size = settings::n_particles;
// int64_t count_size = simulation::work_per_rank;
//
// but that doesn't make sense! The count on a rank is not
// necessarily equal to work_per_rank. For the case of the
// fission source bank, that's true, but for surface source
// creation, there may be a different number of particles
// on each rank.
//
// Hence, I have changed count_size to be the number
// on this rank rather than work_per_rank. Similarly,
// dims_size used to be n_particles, but for surface sources,
// we are not guaranteed to create n_particles at the surface.
// As a result, I've changed this to be the total size of
// the bank across all ranks, which maintains the correct
// behavior for fission bank outputs.
// Set vectors for source bank and starting bank index of each process
vector<int64_t>* bank_index = &simulation::work_index;
vector<SourceSite>* source_bank = &simulation::source_bank;
vector<int64_t> surf_source_index_vector;
vector<SourceSite> surf_source_bank_vector;
// Reset dataspace sizes and vectors for surface source bank
if (surf_source_bank) {
surf_source_index_vector = calculate_surf_source_size();
dims_size = surf_source_index_vector[mpi::n_procs];
count_size = simulation::surf_source_bank.size();
bank_index = &surf_source_index_vector;
// Copy data in a SharedArray into a vector.
surf_source_bank_vector.resize(count_size);
surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
simulation::surf_source_bank.data() + count_size);
source_bank = &surf_source_bank_vector;
}
int64_t dims_size = bank_index.back();
int64_t count_size = bank_index[mpi::rank + 1] - bank_index[mpi::rank];
#ifdef PHDF5
// Set size of total dataspace for all procs and rank
@ -642,7 +631,7 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
hid_t memspace = H5Screate_simple(1, count, nullptr);
// Select hyperslab for this dataspace
hsize_t start[] {static_cast<hsize_t>((*bank_index)[mpi::rank])};
hsize_t start[] {static_cast<hsize_t>(bank_index[mpi::rank])};
H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Set up the property list for parallel writing
@ -650,7 +639,7 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
// Write data to file in parallel
H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank->data());
H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank.data());
// Free resources
H5Sclose(dspace);
@ -669,31 +658,30 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
// Save source bank sites since the array is overwritten below
#ifdef OPENMC_MPI
vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
vector<SourceSite> temp_source {source_bank.begin(), source_bank.end()};
#endif
for (int i = 0; i < mpi::n_procs; ++i) {
// Create memory space
hsize_t count[] {
static_cast<hsize_t>((*bank_index)[i + 1] - (*bank_index)[i])};
hsize_t count[] {static_cast<hsize_t>(bank_index[i + 1] - bank_index[i])};
hid_t memspace = H5Screate_simple(1, count, nullptr);
#ifdef OPENMC_MPI
// Receive source sites from other processes
if (i > 0)
MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
MPI_Recv(source_bank.data(), count[0], mpi::source_site, i, i,
mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// Select hyperslab for this dataspace
dspace = H5Dget_space(dset);
hsize_t start[] {static_cast<hsize_t>((*bank_index)[i])};
hsize_t start[] {static_cast<hsize_t>(bank_index[i])};
H5Sselect_hyperslab(
dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
// Write data to hyperslab
H5Dwrite(
dset, banktype, memspace, dspace, H5P_DEFAULT, (*source_bank).data());
dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank.data());
H5Sclose(memspace);
H5Sclose(dspace);
@ -704,11 +692,11 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
std::copy(temp_source.begin(), temp_source.end(), source_bank.begin());
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
MPI_Send(source_bank.data(), count_size, mpi::source_site, 0, mpi::rank,
mpi::intracomm);
#endif
}

View file

@ -1,5 +1,6 @@
set(TEST_NAMES
test_distribution
test_file_utils
# Add additional unit test files here
)

View file

@ -0,0 +1,28 @@
#include "openmc/file_utils.h"
#include <catch2/catch_test_macros.hpp>
using namespace openmc;
TEST_CASE("Test get_file_extension")
{
REQUIRE(get_file_extension("rememberthealamo.png") == "png");
REQUIRE(get_file_extension("statepoint.20.h5") == "h5");
REQUIRE(get_file_extension("wEiRDNaa_ame.h4") == "h4");
REQUIRE(get_file_extension("has_directory/asdf.20.h5") == "h5");
REQUIRE(get_file_extension("wasssssup_lol") == "");
REQUIRE(get_file_extension("has_directory/secret_file") == "");
}
TEST_CASE("Test dir_exists")
{
// not sure how to test this when running on windows?
REQUIRE(dir_exists("/"));
// if this exists on your system... you deserve for this test to fail
REQUIRE(!dir_exists("/asdfa/asdfasdf/asdgasodgosuihasjkgh/"));
}
TEST_CASE("Test file_exists")
{
// TODO make a file test it exists, delete it
}