diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index 7ae22c0a3..bc5b0a4a9 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -44,6 +44,7 @@ Output Files statepoint source + surface_source summary depletion_results particle_restart diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 588b58b47..53376b04c 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -723,6 +723,42 @@ attributes/sub-elements: *Default*: false +--------------------------- +```` Element +--------------------------- + +The ```` element specifies a surface source file for OpenMC to +read source bank for initializing histories. +This element has the following attributes/sub-elements: + + :path: + Absolute or relative path to a surface source file to read in source bank. + + *Default*: ``surface_source.h5`` in current working directory + +---------------------------- +```` Element +---------------------------- + +The ```` element triggers OpenMC to bank particles crossing +certain surfaces and write out the source bank in a separate file called +``surface_source.h5``. +This element has the following attributes/sub-elements: + + :surface_ids: + A list of integers separated by spaces indicating the unique IDs of surfaces + for which crossing particles will be banked. + + *Default*: None + + :max_particles: + An integer indicating the maximum number of particles to be banked on + specified surfaces per processor. The size of source bank in + ``surface_source.h5`` is limited to this value times the number of + processors. + + *Default*: None + ------------------------------ ```` Element ------------------------------ diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 4cd023e12..86ae8a1e5 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -8,6 +8,10 @@ Normally, source data is stored in a state point file. However, it is possible to request that the source be written separately, in which case the format used is that documented here. +When surface source writing is triggered, a source file named +``surface_source.h5`` is written with only the sources on specified surfaces, +following the same format. + **/** :Attributes: - **filetype** (*char[]*) -- String indicating the type of file. @@ -16,6 +20,7 @@ is that documented here. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, and ``particle``, which represent the - position, direction, energy, weight, delayed group, and particle - type (0=neutron, 1=photon, 2=electron, 3=positron), respectively. + ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, + which represent the position, direction, energy, weight, + delayed group, surface ID, and particle type (0=neutron, 1=photon, + 2=electron, 3=positron), respectively. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 3a2ad1600..92b103ef7 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -52,9 +52,10 @@ The current version of the statepoint file format is 17.0. sum-of-squares for each global tally. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, and ``particle``, which represent the - position, direction, energy, weight, delayed group, and particle - type (0=neutron, 1=photon, 2=electron, 3=positron), respectively. + ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which + represent the position, direction, energy, weight, delayed group, + surface ID, and particle type (0=neutron, 1=photon, 2=electron, + 3=positron), respectively. Only present when `run_mode` is 'eigenvalue'. **/tallies/** diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 7c74125c3..4f6f95c19 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -18,6 +18,8 @@ namespace simulation { extern std::vector source_bank; +extern SharedArray surf_source_bank; + extern SharedArray fission_bank; extern std::vector progeny_per_particle; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index f70faa611..e8002cedb 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -172,6 +172,7 @@ public: double E; double wgt; int delayed_group; + int surf_id; Type particle; int64_t parent_id; int64_t progeny_id; diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 74cbbc132..7030803e9 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -45,6 +45,8 @@ extern "C" bool run_CE; //!< run with continuous-energy data? extern bool source_latest; //!< write latest source at each batch? extern bool source_separate; //!< write source to separate file? extern bool source_write; //!< write source in HDF5 files? +extern bool surf_source_write; //!< write surface source file? +extern bool surf_source_read; //!< read surface source file? extern bool survival_biasing; //!< use survival biasing? extern bool temperature_multipole; //!< use multipole data? extern "C" bool trigger_on; //!< tally triggers enabled? @@ -85,6 +87,8 @@ extern std::vector res_scat_nuclides; //!< Nuclides using res. ups extern RunMode run_mode; //!< Run mode (eigenvalue, fixed src, etc.) extern std::unordered_set sourcepoint_batch; //!< Batches when source should be written extern std::unordered_set statepoint_batch; //!< Batches when state should be written +extern std::unordered_set source_write_surf_id; //!< Surface ids where sources will be written +extern int64_t max_particles; //!< maximum number of particles to be banked on surfaces per process extern TemperatureMethod temperature_method; //!< method for choosing temperatures extern double temperature_tolerance; //!< Tolerance in [K] on choosing temperatures extern double temperature_default; //!< Default T in [K] diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index c50ca355b..53f072b85 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -12,8 +12,9 @@ namespace openmc { void load_state_point(); -void write_source_point(const char* filename); -void write_source_bank(hid_t group_id); +std::vector 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); void read_source_bank(hid_t group_id, std::vector& sites, bool distribute); void write_tally_results_nr(hid_t file_id); void restart_set_keff(); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 1e16c4919..2a62d81c8 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -88,6 +88,7 @@ public: int id_; //!< Unique ID std::string name_; //!< User-defined name std::shared_ptr bc_ {nullptr}; //!< Boundary condition + bool surf_source_ {false}; //!< Activate source banking for the surface? explicit Surface(pugi::xml_node surf_node); Surface(); diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 7be94e03d..6d9a680e3 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -17,6 +17,7 @@ class _Bank(Structure): ('E', c_double), ('wgt', c_double), ('delayed_group', c_int), + ('surf_id', c_int), ('particle', c_int), ('parent_id', c_int64), ('progeny_id', c_int64)] diff --git a/openmc/settings.py b/openmc/settings.py index 4825432a8..dc6a799f5 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -143,6 +143,17 @@ class Settings: Options for writing state points. Acceptable keys are: :batches: list of batches at which to write source + surf_source_read : dict + Options for reading surface source points. Acceptable keys are: + + :path: Path to surface source file (str). + surf_source_write : dict + Options for writing surface source points. Acceptable keys are: + + :surface_ids: List of surface ids at which crossing particles are to be + banked (int) + :max_particles: Maximum number of particles to be banked on surfaces + per process (int) survival_biasing : bool Indicate whether survival biasing is to be used tabular_legendre : dict @@ -229,6 +240,9 @@ class Settings: self._statepoint = {} self._sourcepoint = {} + self._surf_source_read = {} + self._surf_source_write = {} + self._no_reduce = None self._verbosity = None @@ -356,6 +370,14 @@ class Settings: def statepoint(self): return self._statepoint + @property + def surf_source_read(self): + return self._surf_source_read + + @property + def surf_source_write(self): + return self._surf_source_write + @property def no_reduce(self): return self._no_reduce @@ -567,6 +589,35 @@ class Settings: "statepoint options.".format(key)) self._statepoint = statepoint + @surf_source_read.setter + def surf_source_read(self, surf_source_read): + cv.check_type('surface source reading options', surf_source_read, Mapping) + for key, value in surf_source_read.items(): + cv.check_value('surface source reading key', key, + ('path')) + if key == 'path': + cv.check_type('path to surface source file', value, str) + self._surf_source_read = surf_source_read + + @surf_source_write.setter + def surf_source_write(self, surf_source_write): + cv.check_type('surface source writing options', surf_source_write, Mapping) + for key, value in surf_source_write.items(): + cv.check_value('surface source writing key', key, + ('surface_ids', 'max_particles')) + if key == 'surface_ids': + cv.check_type('surface ids for source banking', value, + Iterable, Integral) + for surf_id in value: + cv.check_greater_than('surface id for source banking', + surf_id, 0) + elif key == 'max_particles': + cv.check_type('maximum particle banks on surfaces per process', + value, Integral) + cv.check_greater_than('maximum particle banks on surfaces per process', + value, 0) + self._surf_source_write = surf_source_write + @confidence_intervals.setter def confidence_intervals(self, confidence_intervals): cv.check_type('confidence interval', confidence_intervals, bool) @@ -889,6 +940,24 @@ class Settings: subelement = ET.SubElement(element, "overwrite_latest") subelement.text = str(self._sourcepoint['overwrite']).lower() + def _create_surf_source_read_subelement(self, root): + if self._surf_source_read: + element = ET.SubElement(root, "surf_source_read") + if 'path' in self._surf_source_read: + subelement = ET.SubElement(element, "path") + subelement.text = self._surf_source_read['path'] + + def _create_surf_source_write_subelement(self, root): + if self._surf_source_write: + element = ET.SubElement(root, "surf_source_write") + if 'surface_ids' in self._surf_source_write: + subelement = ET.SubElement(element, "surface_ids") + subelement.text = ' '.join( + str(x) for x in self._surf_source_write['surface_ids']) + if 'max_particles' in self._surf_source_write: + subelement = ET.SubElement(element, "max_particles") + subelement.text = str(self._surf_source_write['max_particles']) + def _create_confidence_intervals(self, root): if self._confidence_intervals is not None: element = ET.SubElement(root, "confidence_intervals") @@ -1151,6 +1220,25 @@ class Settings: value = [int(x) for x in value.split()] self.sourcepoint[key] = value + def _surf_source_read_from_xml_element(self, root): + elem = root.find('surf_source_read') + if elem is not None: + value = get_text(elem, 'path') + if value is not None: + self.surf_source_read['path'] = value + + def _surf_source_write_from_xml_element(self, root): + elem = root.find('surf_source_write') + if elem is not None: + for key in ('surface_ids', 'max_particles'): + value = get_text(elem, key) + if value is not None: + if key == 'surface_ids': + value = [int(x) for x in value.split()] + elif key in ('max_particles'): + value = int(value) + self.surf_source_write[key] = value + def _confidence_intervals_from_xml_element(self, root): text = get_text(root, 'confidence_intervals') if text is not None: @@ -1349,6 +1437,8 @@ class Settings: self._create_output_subelement(root_element) self._create_statepoint_subelement(root_element) self._create_sourcepoint_subelement(root_element) + self._create_surf_source_read_subelement(root_element) + self._create_surf_source_write_subelement(root_element) self._create_confidence_intervals(root_element) self._create_electron_treatment_subelement(root_element) self._create_energy_mode_subelement(root_element) @@ -1422,6 +1512,8 @@ class Settings: settings._output_from_xml_element(root) settings._statepoint_from_xml_element(root) settings._sourcepoint_from_xml_element(root) + settings._surf_source_read_from_xml_element(root) + settings._surf_source_write_from_xml_element(root) settings._confidence_intervals_from_xml_element(root) settings._electron_treatment_from_xml_element(root) settings._energy_mode_from_xml_element(root) diff --git a/openmc/source.py b/openmc/source.py index b5b53167c..3d544f3a1 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -257,17 +257,20 @@ class SourceParticle: Weight of the particle delayed_group : int Delayed group particle was created in (neutrons only) + surf_id : int + Surface ID where particle is at, if any. particle : ParticleType Type of the particle """ def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, wgt=1.0, - delayed_group=0, particle=ParticleType.NEUTRON): + delayed_group=0, surf_id=0, particle=ParticleType.NEUTRON): self.r = tuple(r) self.u = tuple(u) self.E = float(E) self.wgt = float(wgt) self.delayed_group = delayed_group + self.surf_id = surf_id self.particle = particle def to_tuple(self): @@ -280,7 +283,7 @@ class SourceParticle: """ return (self.r, self.u, self.E, self.wgt, - self.delayed_group, self.particle.value) + self.delayed_group, self.surf_id, self.particle.value) def write_source_file(source_particles, filename, **kwargs): @@ -304,6 +307,7 @@ def write_source_file(source_particles, filename, **kwargs): ('E', ' source_bank; +SharedArray surf_source_bank; + // The fission bank is allocated as a SharedArray, rather than a vector, as it will // be shared by all threads in the simulation. It will be allocated to a fixed // maximum capacity in the init_fission_bank() function. Then, Elements will be @@ -37,6 +39,7 @@ std::vector progeny_per_particle; void free_memory_bank() { simulation::source_bank.clear(); + simulation::surf_source_bank.clear(); simulation::fission_bank.clear(); simulation::progeny_per_particle.clear(); } diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 362f9a6e4..58d83f54c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -2,6 +2,7 @@ #include "openmc/cell.h" #include "openmc/constants.h" +#include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/geometry.h" @@ -296,6 +297,10 @@ void load_dagmc_geometry() s->id_ = model::DAG->id_by_index(2, s->dag_index_); s->dagmc_ptr_ = model::DAG; + if (contains(settings::source_write_surf_id, s->id_)) { + s->surf_source_ = true; + } + // set BCs std::string bc_value = DMD.get_surface_property("boundary", surf_handle); to_lower(bc_value); diff --git a/src/initialize.cpp b/src/initialize.cpp index 9581cce0f..6594417de 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -102,22 +102,23 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype Particle::Bank b; - MPI_Aint disp[8]; + MPI_Aint disp[9]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); MPI_Get_address(&b.wgt, &disp[3]); MPI_Get_address(&b.delayed_group, &disp[4]); - MPI_Get_address(&b.particle, &disp[5]); - MPI_Get_address(&b.parent_id, &disp[6]); - MPI_Get_address(&b.progeny_id, &disp[7]); - for (int i = 7; i >= 0; --i) { + MPI_Get_address(&b.surf_id, &disp[5]); + MPI_Get_address(&b.particle, &disp[6]); + MPI_Get_address(&b.parent_id, &disp[7]); + MPI_Get_address(&b.progeny_id, &disp[8]); + for (int i = 8; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(8, blocks, disp, types, &mpi::bank); + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1}; + MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; + MPI_Type_create_struct(9, blocks, disp, types, &mpi::bank); MPI_Type_commit(&mpi::bank); } #endif // OPENMC_MPI diff --git a/src/particle.cpp b/src/particle.cpp index 4453701be..d03511e67 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -413,6 +413,20 @@ Particle::cross_surface() write_message(1, " Crossing surface {}", surf->id_); } + if (surf->surf_source_ && simulation::current_batch == settings::n_batches) { + Particle::Bank site; + site.r = this->r(); + site.u = this->u(); + site.E = this->E_; + site.wgt = this->wgt_; + site.delayed_group = this->delayed_group_; + site.surf_id = surf->id_; + site.particle = this->type_; + site.parent_id = this->id_; + site.progeny_id = this->n_progeny_; + int64_t idx = simulation::surf_source_bank.thread_safe_append(site); + } + // Handle any applicable boundary conditions. if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) { surf->bc_->handle_particle(*this, *surf); diff --git a/src/physics.cpp b/src/physics.cpp index 2636207fb..6e4960513 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -190,6 +190,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) site.wgt = 1. / weight; site.parent_id = p.id_; site.progeny_id = p.n_progeny_++; + site.surf_id = 0; // Sample delayed group and angle/energy for fission reaction sample_fission_neutron(i_nuclide, rx, p.E_, &site, p.current_seed()); diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 94347f864..f9877a3f9 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -139,6 +139,17 @@ element settings { attribute overwrite_latest {xsd:boolean})? }? & + element surf_source_read { + (element path { xsd:string } | attribute path { xsd:string }) + }? & + + element surf_source_write { + (element surface_ids { list { xsd:positiveInteger+ } } | + attribute surface_ids { list { xsd:positiveInteger+ } }) & + (element max_particles { xsd:positiveInteger } | + attribute max_particles { xsd:positiveInteger }) + }? & + element survival_biasing { xsd:boolean }? & element temperature_default { xsd:double }? & diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index ab6dcc328..07b8a6d1d 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -638,6 +638,48 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/settings.cpp b/src/settings.cpp index 5a868d445..9df3b8c46 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -60,6 +60,8 @@ bool run_CE {true}; bool source_latest {false}; bool source_separate {false}; bool source_write {true}; +bool surf_source_write {false}; +bool surf_source_read {false}; bool survival_biasing {false}; bool temperature_multipole {false}; bool trigger_on {false}; @@ -98,6 +100,8 @@ std::vector res_scat_nuclides; RunMode run_mode {RunMode::UNSET}; std::unordered_set sourcepoint_batch; std::unordered_set statepoint_batch; +std::unordered_set source_write_surf_id; +int64_t max_particles; TemperatureMethod temperature_method {TemperatureMethod::NEAREST}; double temperature_tolerance {10.0}; double temperature_default {293.6}; @@ -437,6 +441,20 @@ void read_settings_xml() } } + // Check if the user has specified to read surface source + if (check_for_node(root, "surf_source_read")) { + surf_source_read = true; + // Get surface source read node + xml_node node_ssr = root.child("surf_source_read"); + + std::string path = "surface_source.h5"; + // Check if the user has specified different file for surface source reading + if (check_for_node(node_ssr, "path")) { + path = get_node_value(node_ssr, "path", false, true); + } + model::external_sources.push_back(std::make_unique(path)); + } + // If no source specified, default to isotropic point source at origin with Watt spectrum if (model::external_sources.empty()) { model::external_sources.push_back(std::make_unique( @@ -622,6 +640,26 @@ void read_settings_xml() sourcepoint_batch = statepoint_batch; } + // Check if the user has specified to write surface source + if (check_for_node(root, "surf_source_write")) { + surf_source_write = true; + // Get surface source write node + xml_node node_ssw = root.child("surf_source_write"); + + // Determine surface ids at which crossing particles are to be banked + if (check_for_node(node_ssw, "surface_ids")) { + auto temp = get_node_array(node_ssw, "surface_ids"); + for (const auto& b : temp) { + source_write_surf_id.insert(b); + } + } + + // Get maximum number of particles to be banked per surface + if (check_for_node(node_ssw, "max_particles")) { + max_particles = std::stoi(get_node_value(node_ssw, "max_particles")); + } + } + // If source is not seperate and is to be written out in the statepoint file, // make sure that the sourcepoint batch numbers are contained in the // statepoint list @@ -794,6 +832,7 @@ void read_settings_xml() void free_memory_settings() { settings::statepoint_batch.clear(); settings::sourcepoint_batch.clear(); + settings::source_write_surf_id.clear(); settings::res_scat_nuclides.clear(); } diff --git a/src/simulation.cpp b/src/simulation.cpp index 1e35adeee..303575130 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -79,10 +79,8 @@ int openmc_simulation_init() // Determine how much work each process should do calculate_work(); - // Allocate source and fission banks for eigenvalue simulations - if (settings::run_mode == RunMode::EIGENVALUE) { - allocate_banks(); - } + // Allocate source, fission and surface source banks. + allocate_banks(); // If doing an event-based simulation, intialize the particle buffer // and event queues @@ -281,10 +279,19 @@ std::vector work_index; void allocate_banks() { - // Allocate source bank - simulation::source_bank.resize(simulation::work_per_rank); - // Allocate fission bank - init_fission_bank(3*simulation::work_per_rank); + if (settings::run_mode == RunMode::EIGENVALUE) { + // Allocate source bank + simulation::source_bank.resize(simulation::work_per_rank); + + // Allocate fission bank + init_fission_bank(3*simulation::work_per_rank); + } + + if (settings::surf_source_write) { + // Allocate surface source bank + simulation::surf_source_bank.reserve(settings::max_particles); + } + } void initialize_batch() @@ -375,6 +382,12 @@ void finalize_batch() write_source_point(filename.c_str()); } } + + // Write out surface source if requested. + if (settings::surf_source_write && simulation::current_batch == settings::n_batches) { + auto filename = settings::path_output + "surface_source.h5"; + write_source_point(filename.c_str(), true); + } } void initialize_generation() @@ -439,7 +452,10 @@ void finalize_generation() void initialize_history(Particle& p, int64_t index_source) { // set defaults - if (settings::run_mode == RunMode::FIXED_SOURCE) { + if (settings::run_mode == RunMode::EIGENVALUE) { + // set defaults for eigenvalue simulations from primary bank + p.from_source(&simulation::source_bank[index_source - 1]); + } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + simulation::work_index[mpi::rank] + index_source; @@ -447,9 +463,6 @@ void initialize_history(Particle& p, int64_t index_source) // sample from external source distribution or custom library then set auto site = sample_external_source(&seed); p.from_source(&site); - } else if (settings::run_mode == RunMode::EIGENVALUE) { - // set defaults for eigenvalue simulations from primary bank - p.from_source(&simulation::source_bank[index_source - 1]); } p.current_work_ = index_source; diff --git a/src/source.cpp b/src/source.cpp index 87644f054..2bbe7c321 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -225,6 +225,8 @@ Particle::Bank IndependentSource::sample(uint64_t* seed) const // Set delayed group site.delayed_group = 0; + // Set surface ID + site.surf_id = 0; return site; } @@ -342,7 +344,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); + write_source_bank(file_id, false); file_close(file_id); } } diff --git a/src/state_point.cpp b/src/state_point.cpp index 66706a676..48d87d403 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -311,7 +311,7 @@ openmc_statepoint_write(const char* filename, bool* write_source) // Write the source bank if desired if (write_source_) { if (mpi::master || parallel) file_id = file_open(filename_, 'a', true); - write_source_bank(file_id); + write_source_bank(file_id, false); if (mpi::master || parallel) file_close(file_id); } @@ -518,14 +518,38 @@ hid_t h5banktype() { H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE); H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE); H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT); + H5Tinsert(banktype, "surf_id", HOFFSET(Particle::Bank, surf_id), H5T_NATIVE_INT); H5Tinsert(banktype, "particle", HOFFSET(Particle::Bank, particle), H5T_NATIVE_INT); H5Tclose(postype); return banktype; } +std::vector calculate_surf_source_size() +{ + std::vector surf_source_index; + surf_source_index.reserve(mpi::n_procs + 1); + +#ifdef OPENMC_MPI + surf_source_index.resize(mpi::n_procs); + std::vector 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) +write_source_point(const char* filename, bool surf_source_bank) { // 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 @@ -555,29 +579,54 @@ write_source_point(const char* filename) } // Get pointer to source bank and write to file - write_source_bank(file_id); + write_source_bank(file_id, surf_source_bank); if (mpi::master || parallel) file_close(file_id); } void -write_source_bank(hid_t group_id) +write_source_bank(hid_t group_id, bool surf_source_bank) { 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; + + // Set vectors for source bank and starting bank index of each process + std::vector* bank_index = &simulation::work_index; + std::vector* source_bank = &simulation::source_bank; + std::vector surf_source_index_vector; + std::vector 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; + } + #ifdef PHDF5 // Set size of total dataspace for all procs and rank - hsize_t dims[] {static_cast(settings::n_particles)}; + hsize_t dims[] {static_cast(dims_size)}; hid_t dspace = H5Screate_simple(1, dims, nullptr); hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); // Create another data space but for each proc individually - hsize_t count[] {static_cast(simulation::work_per_rank)}; + hsize_t count[] {static_cast(count_size)}; hid_t memspace = H5Screate_simple(1, count, nullptr); // Select hyperslab for this dataspace - hsize_t start[] {static_cast(simulation::work_index[mpi::rank])}; + hsize_t start[] {static_cast((*bank_index)[mpi::rank])}; H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); // Set up the property list for parallel writing @@ -585,7 +634,7 @@ write_source_bank(hid_t group_id) H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE); // Write data to file in parallel - H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data()); + H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank->data()); // Free resources H5Sclose(dspace); @@ -597,38 +646,35 @@ write_source_bank(hid_t group_id) if (mpi::master) { // Create dataset big enough to hold all source sites - hsize_t dims[] {static_cast(settings::n_particles)}; + hsize_t dims[] {static_cast(dims_size)}; hid_t dspace = H5Screate_simple(1, dims, nullptr); hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - // Save source bank sites since the souce_bank array is overwritten below + // Save source bank sites since the array is overwritten below #ifdef OPENMC_MPI - std::vector temp_source {simulation::source_bank.begin(), - simulation::source_bank.begin() + simulation::work_per_rank}; + std::vector 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(simulation::work_index[i+1] - - simulation::work_index[i])}; + hsize_t count[] {static_cast((*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(simulation::source_bank.data(), count[0], mpi::bank, i, i, + MPI_Recv(source_bank->data(), count[0], mpi::bank, i, i, mpi::intracomm, MPI_STATUS_IGNORE); #endif // Select hyperslab for this dataspace dspace = H5Dget_space(dset); - hsize_t start[] {static_cast(simulation::work_index[i])}; + hsize_t start[] {static_cast((*bank_index)[i])}; H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr); // Write data to hyperslab - H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, - simulation::source_bank.data()); + H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, (*source_bank).data()); H5Sclose(memspace); H5Sclose(dspace); @@ -639,11 +685,11 @@ write_source_bank(hid_t group_id) #ifdef OPENMC_MPI // Restore state of source bank - std::copy(temp_source.begin(), temp_source.end(), simulation::source_bank.begin()); + std::copy(temp_source.begin(), temp_source.end(), source_bank->begin()); #endif } else { #ifdef OPENMC_MPI - MPI_Send(simulation::source_bank.data(), simulation::work_per_rank, mpi::bank, + MPI_Send(source_bank->data(), count_size, mpi::bank, 0, mpi::rank, mpi::intracomm); #endif } diff --git a/src/surface.cpp b/src/surface.cpp index e0f2787b6..81d96d83a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -8,6 +8,7 @@ #include #include +#include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/dagmc.h" #include "openmc/hdf5_interface.h" @@ -115,6 +116,9 @@ Surface::Surface(pugi::xml_node surf_node) { if (check_for_node(surf_node, "id")) { id_ = std::stoi(get_node_value(surf_node, "id")); + if (contains(settings::source_write_surf_id, id_)) { + surf_source_ = true; + } } else { fatal_error("Must specify id of surface in geometry XML file."); } diff --git a/tests/regression_tests/surface_source/__init__.py b/tests/regression_tests/surface_source/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/surface_source/inputs_true_read.dat b/tests/regression_tests/surface_source/inputs_true_read.dat new file mode 100644 index 000000000..3a122c737 --- /dev/null +++ b/tests/regression_tests/surface_source/inputs_true_read.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + surface_source_true.h5 + + 1 + + + + + 3 + + + 1 + flux + + diff --git a/tests/regression_tests/surface_source/inputs_true_write.dat b/tests/regression_tests/surface_source/inputs_true_write.dat new file mode 100644 index 000000000..48dde9900 --- /dev/null +++ b/tests/regression_tests/surface_source/inputs_true_write.dat @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + 0 0 0 + + + + 1 + 1000 + + 1 + + + + + 3 + + + 1 + flux + + diff --git a/tests/regression_tests/surface_source/results_true.dat b/tests/regression_tests/surface_source/results_true.dat new file mode 100644 index 000000000..558217673 --- /dev/null +++ b/tests/regression_tests/surface_source/results_true.dat @@ -0,0 +1,3 @@ +tally 1: +5.000000E+00 +2.500000E+00 diff --git a/tests/regression_tests/surface_source/surface_source_true.h5 b/tests/regression_tests/surface_source/surface_source_true.h5 new file mode 100644 index 000000000..8aefd8159 Binary files /dev/null and b/tests/regression_tests/surface_source/surface_source_true.h5 differ diff --git a/tests/regression_tests/surface_source/test.py b/tests/regression_tests/surface_source/test.py new file mode 100644 index 000000000..492ef192e --- /dev/null +++ b/tests/regression_tests/surface_source/test.py @@ -0,0 +1,125 @@ +import os + +import h5py +import numpy as np +import pytest +import openmc + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(request): + openmc.reset_auto_ids() + marker = request.node.get_closest_marker("surf_source_op") + surf_source_op = marker.args[0] + + openmc_model = openmc.model.Model() + + # Materials + # None + + # Geometry + # Concentric void spheres + # - Innermost sphere to bank surface sources + # - Second shell to tally cell flux + # - Outermost sphere as vacuum boundary + sph_1 = openmc.Sphere(r=1.0) # Surface to bank/write sources. + sph_2 = openmc.Sphere(r=2.0) + sph_3 = openmc.Sphere(r=2.5) + sph_4 = openmc.Sphere(r=4.0, boundary_type='vacuum') + cell_1 = openmc.Cell(region=-sph_1) + cell_2 = openmc.Cell(region=+sph_1&-sph_2) + cell_3 = openmc.Cell(region=+sph_2&-sph_3) # Cell to tally flux. + cell_4 = openmc.Cell(region=+sph_3&-sph_4) + root = openmc.Universe(cells=[cell_1, cell_2, cell_3, cell_4]) + openmc_model.geometry = openmc.Geometry(root) + + # Settings + openmc_model.settings.run_mode = 'fixed source' + openmc_model.settings.particles = 1000 + openmc_model.settings.batches = 10 + openmc_model.settings.seed = 1 + + if surf_source_op == 'write': + point = openmc.stats.Point((0, 0, 0)) + pt_src = openmc.Source(space=point) + openmc_model.settings.source = pt_src + + openmc_model.settings.surf_source_write = {'surface_ids': [1], + 'max_particles': 1000} + elif surf_source_op == 'read': + openmc_model.settings.surf_source_read = {'path': 'surface_source_true.h5'} + + # Tallies + tal = openmc.Tally() + cell_filter = openmc.CellFilter(cell_3) + tal.filters = [cell_filter] + tal.scores = ['flux'] + openmc_model.tallies.append(tal) + + return openmc_model + + +class SurfaceSourceTestHarness(PyAPITestHarness): + def _test_output_created(self): + """Make sure surface_source.h5 has also been created.""" + super()._test_output_created() + # Check if 'surface_source.h5' has been created. + if self._model.settings.surf_source_write: + assert os.path.exists('surface_source.h5'), \ + 'Surface source file does not exist.' + + def _compare_output(self): + """Make sure the current surface_source.h5 agree with the reference.""" + if self._model.settings.surf_source_write: + with h5py.File("surface_source_true.h5", 'r') as f: + source_true = f['source_bank'][()] + # Convert dtye from mixed to a float for comparison assertion + source_true.dtype = 'float64' + with h5py.File("surface_source.h5", 'r') as f: + source_test = f['source_bank'][()] + # Convert dtye from mixed to a float for comparison assertion + source_test.dtype = 'float64' + np.testing.assert_allclose(np.sort(source_true), + np.sort(source_test), + atol=1e-07) + + def execute_test(self): + """Build input XMLs, run OpenMC, check output and results.""" + try: + self._build_inputs() + inputs = self._get_inputs() + self._write_inputs(inputs) + self._compare_inputs() + self._run_openmc() + self._test_output_created() + self._compare_output() + results = self._get_results() + self._write_results(results) + self._compare_results() + finally: + self._cleanup() + + def _cleanup(self): + """Delete statepoints, tally, and test files.""" + super()._cleanup() + fs = 'surface_source.h5' + if os.path.exists(fs): + os.remove(fs) + + +@pytest.mark.surf_source_op('write') +def test_surface_source_write(model): + harness = SurfaceSourceTestHarness('statepoint.10.h5', + model, + 'inputs_true_write.dat') + harness.main() + + +@pytest.mark.surf_source_op('read') +def test_surface_source_read(model): + harness = SurfaceSourceTestHarness('statepoint.10.h5', + model, + 'inputs_true_read.dat') + harness.main() diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index c3ff72879..7a3357d43 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -20,6 +20,8 @@ def test_export_to_xml(run_in_tmpdir): s.sourcepoint = {'batches': [50, 150, 500, 1000], 'separate': True, 'write': True, 'overwrite': True} s.statepoint = {'batches': [50, 150, 500, 1000]} + s.surf_source_read = {'path': 'surface_source_1.h5'} + s.surf_source_write = {'surface_ids': [2], 'max_particles': 200} s.confidence_intervals = True s.ptables = True s.seed = 17 @@ -76,6 +78,8 @@ def test_export_to_xml(run_in_tmpdir): assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True, 'write': True, 'overwrite': True} assert s.statepoint == {'batches': [50, 150, 500, 1000]} + assert s.surf_source_read == {'path': 'surface_source_1.h5'} + assert s.surf_source_write == {'surface_ids': [2], 'max_particles': 200} assert s.confidence_intervals assert s.ptables assert s.seed == 17