diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 93e2be669..303c1e8be 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -7,11 +7,6 @@ #include "openmc/particle.h" #include "openmc/position.h" -// Without an explicit instantiation of vector, the Intel compiler -// will complain about the threadprivate directive on filter_matches. Note that -// this has to happen *outside* of the openmc namespace -extern template class std::vector; - namespace openmc { //============================================================================== @@ -21,13 +16,11 @@ namespace openmc { namespace simulation { extern std::vector source_bank; -extern std::vector fission_bank; -extern std::vector secondary_bank; -#ifdef _OPENMP -extern std::vector master_fission_bank; -#endif -#pragma omp threadprivate(fission_bank, secondary_bank) +extern std::unique_ptr fission_bank; +extern int64_t fission_bank_length; +extern int64_t fission_bank_max; +extern std::vector progeny_per_particle; } // namespace simulation @@ -35,8 +28,12 @@ extern std::vector master_fission_bank; // Non-member functions //============================================================================== +void sort_fission_bank(); + void free_memory_bank(); +void init_fission_bank(int64_t max); + } // namespace openmc #endif // OPENMC_BANK_H diff --git a/include/openmc/cell.h b/include/openmc/cell.h index e5797872e..a5e530d7d 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -113,7 +113,7 @@ public: //! Find the oncoming boundary of this cell. virtual std::pair - distance(Position r, Direction u, int32_t on_surface) const = 0; + distance(Position r, Direction u, int32_t on_surface, Particle* p) const = 0; //! Write all information needed to reconstruct the cell to an HDF5 group. //! \param group_id An HDF5 group id. @@ -204,7 +204,7 @@ public: contains(Position r, Direction u, int32_t on_surface) const; std::pair - distance(Position r, Direction u, int32_t on_surface) const; + distance(Position r, Direction u, int32_t on_surface, Particle* p) const; void to_hdf5(hid_t group_id) const; @@ -248,7 +248,7 @@ public: bool contains(Position r, Direction u, int32_t on_surface) const; std::pair - distance(Position r, Direction u, int32_t on_surface) const; + distance(Position r, Direction u, int32_t on_surface, Particle* p) const; BoundingBox bounding_box() const; diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index d3012e3a0..859a229f0 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -14,14 +14,6 @@ extern "C" const bool dagmc_enabled; namespace openmc { -namespace simulation { - -extern moab::DagMC::RayHistory history; //!< facet history for DagMC particles -extern Direction last_dir; //!< last direction passed to DagMC's ray_fire -#pragma omp threadprivate(history, last_dir) - -} - namespace model { extern moab::DagMC* DAG; } diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index d29f404b9..0e2d2401d 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -42,14 +42,6 @@ void calculate_generation_keff(); //! generations. It also broadcasts the value from the master process. void calculate_average_keff(); -#ifdef _OPENMP -//! Join threadprivate fission banks into a single fission bank -//! -//! Note that this operation is necessarily sequential to preserve the order of -//! the bank when using varying numbers of threads. -void join_bank_from_threads(); -#endif - //! Calculates a minimum variance estimate of k-effective //! //! The minimum variance estimate is based on a linear combination of the diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index e386348c4..ef226db66 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -24,17 +24,6 @@ extern std::vector overlap_check_count; } // namespace model -//============================================================================== -// Information about nearest boundary crossing -//============================================================================== - -struct BoundaryInfo { - double distance {INFINITY}; //!< distance to nearest boundary - int surface_index {0}; //!< if boundary is surface, index in surfaces vector - int coord_level; //!< coordinate level after crossing boundary - std::array lattice_translation {}; //!< which way lattice indices will change -}; - //============================================================================== //! Check two distances by coincidence tolerance //============================================================================== diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 34782b5e3..bac4885ca 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -161,7 +161,7 @@ public: //! \param[in] bank Array of bank sites //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xtensor count_sites(const std::vector& bank, + xt::xtensor count_sites(const Particle::Bank* bank, int64_t length, bool* outside) const; // Data members diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 1b049c2e2..8143ee5a3 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -14,6 +14,12 @@ #include "openmc/constants.h" #include "openmc/position.h" #include "openmc/random_lcg.h" +#include "openmc/tallies/filter_match.h" + +#ifdef DAGMC +#include "DagMC.hpp" +#endif + namespace openmc { @@ -131,6 +137,17 @@ struct MacroXS { double pair_production; //!< macroscopic pair production xs }; +//============================================================================== +// Information about nearest boundary crossing +//============================================================================== + +struct BoundaryInfo { + double distance {INFINITY}; //!< distance to nearest boundary + int surface_index {0}; //!< if boundary is surface, index in surfaces vector + int coord_level; //!< coordinate level after crossing boundary + std::array lattice_translation {}; //!< which way lattice indices will change +}; + //============================================================================ //! State of a particle being transported through geometry //============================================================================ @@ -146,6 +163,12 @@ public: }; //! Saved ("banked") state of a particle + //! NOTE: This structure's MPI type is built in initialize_mpi() of + //! initialize.cpp. Any changes made to the struct here must also be + //! made when building the Bank MPI type in initialize_mpi(). + //! NOTE: This structure is also used on the python side, and is defined + //! in lib/core.py. Changes made to the type here must also be made to the + //! python defintion. struct Bank { Position r; Direction u; @@ -153,6 +176,15 @@ public: double wgt; int delayed_group; Type particle; + int64_t parent_id; + int64_t progeny_id; + }; + + //! Saved ("banked") state of a particle, for nu-fission tallying + struct NuBank { + double E; //!< particle energy + double wgt; //!< particle weight + int delayed_group; //!< particle delayed group }; //========================================================================== @@ -199,8 +231,13 @@ public: //! \param src Source site data void from_source(const Bank* src); - //! Transport a particle from birth to death - void transport(); + // Coarse-grained particle events + void event_calculate_xs(); + void event_advance(); + void event_cross_surface(); + void event_collide(); + void event_revive_from_secondary(); + void event_death(); //! Cross a surface and handle boundary conditions void cross_surface(); @@ -276,10 +313,13 @@ public: //!< sites banked // Indices for various arrays - int surface_ {0}; //!< index for surface particle is on + int surface_ {0}; //!< index for surface particle is on int cell_born_ {-1}; //!< index for cell particle was born in int material_ {-1}; //!< index for current material int material_last_ {-1}; //!< index for last material + + // Boundary information + BoundaryInfo boundary_; // Temperature of current cell double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV @@ -294,6 +334,39 @@ public: // Current PRNG state uint64_t seeds_[N_STREAMS]; // current seeds int stream_; // current RNG stream + + // Secondary particle bank + std::vector secondary_bank_; + + int64_t current_work_; // current work index + + std::vector flux_derivs_; // for derivatives for this particle + + std::vector filter_matches_; // tally filter matches + + std::vector> tracks_; // tracks for outputting to file + + std::vector nu_bank_; // bank of most recently fissioned particles + + // Global tally accumulators + double keff_tally_absorption_ {0.0}; + double keff_tally_collision_ {0.0}; + double keff_tally_tracklength_ {0.0}; + double keff_tally_leakage_ {0.0}; + + bool trace_ {false}; //!< flag to show debug information + + double collision_distance_; // distance to particle's next closest collision + + int n_event_ {0}; // number of events executed in this particle's history + + // DagMC state variables + #ifdef DAGMC + moab::DagMC::RayHistory history_; + Direction last_dir_; + #endif + + int64_t n_progeny_ {0}; // Number of progeny produced by this particle }; } // namespace openmc diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 048ba7998..ac83dbc70 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -48,8 +48,7 @@ int sample_nuclide(Particle* p); //! Determine the average total, prompt, and delayed neutrons produced from //! fission and creates appropriate bank sites. -void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, - std::vector& bank); +void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx); int sample_element(Particle* p); diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index f246cb6ec..5cb7226d9 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -33,9 +33,8 @@ scatter(Particle* p); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. //! \param p Particle to operate on -//! \param bank The particle bank to populate void -create_fission_sites(Particle* p, std::vector& bank); +create_fission_sites(Particle* p); //! \brief Handles an absorption event //! \param p Particle to operate on diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 8efac4c7f..5bfba138c 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -24,7 +24,6 @@ namespace simulation { extern "C" int current_batch; //!< current batch extern "C" int current_gen; //!< current fission generation -extern "C" int64_t current_work; //!< index in source back of current particle extern "C" bool initialized; //!< has simulation been initialized? extern "C" double keff; //!< average k over batches extern "C" double keff_std; //!< standard deviation of average k @@ -46,11 +45,6 @@ extern const RegularMesh* ufs_mesh; extern std::vector k_generation; extern std::vector work_index; -// Threadprivate variables -extern "C" bool trace; //!< flag to show debug information - -#pragma omp threadprivate(current_work, trace) - } // namespace simulation //============================================================================== @@ -69,8 +63,12 @@ void initialize_batch(); //! Initialize a fission generation void initialize_generation(); +//! Full initialization of a particle history void initialize_history(Particle* p, int64_t index_source); +//! Helper function for initialize_history() that is called independently elsewhere +void initialize_history_partial(Particle* p); + //! Finalize a batch //! //! Handles synchronization and accumulation of tallies, calculation of Shannon @@ -90,6 +88,13 @@ void broadcast_results(); void free_memory_simulation(); +//! Simulate a single particle history (and all generated secondary particles, +//! if enabled), from birth to death +void transport_history_based_single_particle(Particle& p); + +//! Simulate all particle histories using history-based parallelism +void transport_history_based(); + } // namespace openmc #endif // OPENMC_SIMULATION_H diff --git a/include/openmc/surface.h b/include/openmc/surface.h index c72fa8ef3..3ba147a5f 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -11,6 +11,7 @@ #include "pugixml.hpp" #include "openmc/constants.h" +#include "openmc/particle.h" #include "openmc/position.h" #include "dagmc.h" @@ -112,8 +113,9 @@ public: //! Determine the direction of a ray reflected from the surface. //! \param[in] r The point at which the ray is incident. //! \param[in] u Incident direction of the ray + //! \param[inout] p Pointer to the particle //! \return Outgoing direction of the ray - virtual Direction reflect(Position r, Direction u) const; + virtual Direction reflect(Position r, Direction u, Particle* p) const; virtual Direction diffuse_reflect(Position r, Direction u, uint64_t* seed) const; @@ -170,7 +172,7 @@ public: double evaluate(Position r) const; double distance(Position r, Direction u, bool coincident) const; Direction normal(Position r) const; - Direction reflect(Position r, Direction u) const; + Direction reflect(Position r, Direction u, Particle* p) const; void to_hdf5(hid_t group_id) const; diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index be998e80b..5cfbfc712 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -27,7 +27,6 @@ struct TallyDerivative { int id; //!< User-defined identifier int diff_material; //!< Material this derivative is applied to int diff_nuclide; //!< Nuclide this material is applied to - double flux_deriv; //!< Derivative of the current particle's weight TallyDerivative() {} explicit TallyDerivative(pugi::xml_node node); @@ -54,16 +53,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, //! further tallies are scored. // //! \param p The particle being tracked -void score_collision_derivative(const Particle* p); +void score_collision_derivative(Particle* p); //! Adjust diff tally flux derivatives for a particle tracking event. // //! \param p The particle being tracked //! \param distance The distance in [cm] traveled by the particle -void score_track_derivative(const Particle* p, double distance); - -//! Set the flux derivatives on differential tallies to zero. -void zero_flux_derivs(); +void score_track_derivative(Particle* p, double distance); } // namespace openmc @@ -71,15 +67,10 @@ void zero_flux_derivs(); // Global variables //============================================================================== -// Explicit vector template specialization declaration of threadprivate variable -// outside of the openmc namespace for the picky Intel compiler. -extern template class std::vector; - namespace openmc { namespace model { extern std::vector tally_derivs; -#pragma omp threadprivate(tally_derivs) extern std::unordered_map tally_deriv_map; } // namespace model diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index b95d95bbb..3fd0e97f3 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -12,31 +12,10 @@ #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" +#include "openmc/tallies/filter_match.h" #include "pugixml.hpp" -namespace openmc { - -//============================================================================== -//! Stores bins and weights for filtered tally events. -//============================================================================== - -class FilterMatch -{ -public: - std::vector bins_; - std::vector weights_; - int i_bin_; - bool bins_present_ {false}; -}; - -} // namespace openmc - -// Without an explicit instantiation of vector, the Intel compiler -// will complain about the threadprivate directive on filter_matches. Note that -// this has to happen *outside* of the openmc namespace -extern template class std::vector; - namespace openmc { //============================================================================== @@ -127,13 +106,6 @@ private: // Global variables //============================================================================== -namespace simulation { - -extern std::vector filter_matches; -#pragma omp threadprivate(filter_matches) - -} // namespace simulation - namespace model { extern "C" int32_t n_filters; extern std::vector> tally_filters; diff --git a/include/openmc/tallies/filter_match.h b/include/openmc/tallies/filter_match.h new file mode 100644 index 000000000..a4ab2e407 --- /dev/null +++ b/include/openmc/tallies/filter_match.h @@ -0,0 +1,21 @@ +#ifndef OPENMC_TALLIES_FILTERMATCH_H +#define OPENMC_TALLIES_FILTERMATCH_H + + +namespace openmc { + +//============================================================================== +//! Stores bins and weights for filtered tally events. +//============================================================================== + +class FilterMatch +{ +public: + std::vector bins_; + std::vector weights_; + int i_bin_; + bool bins_present_ {false}; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTERMATCH_H diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index e8f2ff5e5..96dc07fb9 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -160,17 +160,10 @@ namespace simulation { extern "C" int32_t n_realizations; } -// It is possible to protect accumulate operations on global tallies by using an -// atomic update. However, when multiple threads accumulate to the same global -// tally, it can cause a higher cache miss rate due to invalidation. Thus, we -// use threadprivate variables to accumulate global tallies and then reduce at -// the end of a generation. extern double global_tally_absorption; extern double global_tally_collision; extern double global_tally_tracklength; extern double global_tally_leakage; -#pragma omp threadprivate(global_tally_absorption, global_tally_collision, \ - global_tally_tracklength, global_tally_leakage) //============================================================================== // Non-member functions diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 2c5a29c17..224a9ef14 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -2,6 +2,7 @@ #define OPENMC_TALLIES_TALLY_SCORING_H #include "openmc/particle.h" +#include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" namespace openmc { @@ -22,12 +23,13 @@ class FilterBinIter public: //! Construct an iterator over bins that match a given particle's state. - FilterBinIter(const Tally& tally, const Particle* p); + FilterBinIter(const Tally& tally, Particle* p); //! Construct an iterator over all filter bin combinations. // //! \param end if true, the returned iterator indicates the end of a loop. - FilterBinIter(const Tally& tally, bool end); + FilterBinIter(const Tally& tally, bool end, + std::vector* particle_filter_matches); bool operator==(const FilterBinIter& other) const {return index_ == other.index_;} @@ -39,6 +41,8 @@ public: int index_ {1}; double weight_ {1.}; + + std::vector& filter_matches_; private: void compute_index_weight(); @@ -73,7 +77,7 @@ void score_analog_tally_ce(Particle* p); //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_mg(const Particle* p); +void score_analog_tally_mg(Particle* p); //! Score tallies using a tracklength estimate of the flux. // @@ -89,7 +93,7 @@ void score_tracklength_tally(Particle* p, double distance); // //! \param p The particle being tracked //! \param tallies A vector of tallies to score to -void score_surface_tally(const Particle* p, const std::vector& tallies); +void score_surface_tally(Particle* p, const std::vector& tallies); } // namespace openmc diff --git a/include/openmc/timer.h b/include/openmc/timer.h index 9149c4ef7..b2d5f0005 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -21,6 +21,7 @@ extern Timer time_finalize; extern Timer time_inactive; extern Timer time_initialize; extern Timer time_read_xs; +extern Timer time_sample_source; extern Timer time_tallies; extern Timer time_total; extern Timer time_transport; diff --git a/include/openmc/track_output.h b/include/openmc/track_output.h index f32529bc2..fc65a2eaa 100644 --- a/include/openmc/track_output.h +++ b/include/openmc/track_output.h @@ -9,9 +9,9 @@ namespace openmc { // Non-member functions //============================================================================== -void add_particle_track(); -void write_particle_track(const Particle& p); -void finalize_particle_track(const Particle& p); +void add_particle_track(Particle& p); +void write_particle_track(Particle& p); +void finalize_particle_track(Particle& p); } // namespace openmc diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 83e773f20..de2e5e924 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1073,11 +1073,15 @@ class Reaction(EqualityMixin): if i_reaction < ace.nxs[5] + 1: # Check if angular distribution data exist loc = int(ace.xss[ace.jxs[8] + i_reaction]) - if loc <= 0: - # Angular distribution is either given as part of a product - # angle-energy distribution or is not given at all (in which - # case isotropic scattering is assumed) + if loc < 0: + # Angular distribution is given as part of a product + # angle-energy distribution angle_dist = None + elif loc == 0: + # Angular distribution is isotropic + energy = [0.0, grid[-1]] + mu = Uniform(-1., 1.) + angle_dist = AngleDistribution(energy, [mu, mu]) else: angle_dist = AngleDistribution.from_ace(ace, ace.jxs[9], loc) diff --git a/openmc/geometry.py b/openmc/geometry.py index 5424f30e5..b911b0165 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -77,15 +77,22 @@ class Geometry(object): if universe.id in volume_calc.volumes: universe.add_volume_information(volume_calc) - def export_to_xml(self, path='geometry.xml'): + def export_to_xml(self, path='geometry.xml', remove_surfs=False): """Export geometry to an XML file. Parameters ---------- path : str Path to file to write. Defaults to 'geometry.xml'. + remove_surfs : bool + Whether or not to remove redundant surfaces from the geometry when + exporting """ + # Find and remove redundant surfaces from the geometry + if remove_surfs: + self.remove_redundant_surfaces() + # Create XML representation root_element = ET.Element("geometry") self.root_universe.create_xml_subelement(root_element, memo=set()) @@ -382,6 +389,26 @@ class Geometry(object): surfaces = cell.region.get_surfaces(surfaces) return surfaces + def get_redundant_surfaces(self): + """Return all of the topologically redundant surface ids + + Returns + ------- + dict + Dictionary whose keys are the ID of a redundant surface and whose + values are the topologically equivalent :class:`openmc.Surface` + that should replace it. + + """ + tally = defaultdict(list) + for surf in self.get_all_surfaces().values(): + coeffs = tuple(surf._coefficients[k] for k in surf._coeff_keys) + key = (surf._type,) + coeffs + tally[key].append(surf) + return {replace.id: keep + for keep, *redundant in tally.values() + for replace in redundant} + def get_materials_by_name(self, name, case_sensitive=False, matching=False): """Return a list of materials with matching names. @@ -579,6 +606,17 @@ class Geometry(object): return sorted(lattices, key=lambda x: x.id) + def remove_redundant_surfaces(self): + """Remove redundant surfaces from the geometry""" + + # Get redundant surfaces + redundant_surfaces = self.get_redundant_surfaces() + + # Iterate through all cells contained in the geometry + for cell in self.get_all_cells().values(): + # Recursively remove redundant surfaces from regions + cell.region.remove_redundant_surfaces(redundant_surfaces) + def determine_paths(self, instances_only=False): """Determine paths through CSG tree for cells and materials. diff --git a/openmc/lib/core.py b/openmc/lib/core.py index cba96fa4b..dcd896fa6 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -20,7 +20,9 @@ class _Bank(Structure): ('E', c_double), ('wgt', c_double), ('delayed_group', c_int), - ('particle', c_int)] + ('particle', c_int), + ('parent_id', c_int64), + ('progeny_id', c_int64)] # Define input type for numpy arrays that will be passed into C++ functions diff --git a/openmc/region.py b/openmc/region.py index a751e1ebf..2b2fce14d 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -45,8 +45,7 @@ class Region(metaclass=ABCMeta): return not self == other def get_surfaces(self, surfaces=None): - """ - Recursively find all the surfaces referenced by a region and return them + """Recursively find all surfaces referenced by a region and return them Parameters ---------- @@ -65,6 +64,19 @@ class Region(metaclass=ABCMeta): surfaces = region.get_surfaces(surfaces) return surfaces + def remove_redundant_surfaces(self, redundant_surfaces): + """Recursively remove all redundant surfaces referenced by this region + + Parameters + ---------- + redundant_surfaces : dict + Dictionary mapping redundant surface IDs to class:`openmc.Surface` + instances that should replace them. + + """ + for region in self: + region.remove_redundant_surfaces(redundant_surfaces) + @staticmethod def from_expression(expression, surfaces): """Generate a region given an infix expression. @@ -597,8 +609,7 @@ class Complement(Region): return temp_region.bounding_box def get_surfaces(self, surfaces=None): - """ - Recursively find and return all the surfaces referenced by the node + """Recursively find and return all the surfaces referenced by the node Parameters ---------- @@ -617,6 +628,19 @@ class Complement(Region): surfaces = region.get_surfaces(surfaces) return surfaces + def remove_redundant_surfaces(self, redundant_surfaces): + """Recursively remove all redundant surfaces referenced by this region + + Parameters + ---------- + redundant_surfaces : dict + Dictionary mapping redundant surface IDs to class:`openmc.Surface` + instances that should replace them. + + """ + for region in self.node: + region.remove_redundant_surfaces(redundant_surfaces) + def clone(self, memo=None): """Create a copy of this region - each of the surfaces in the complement's node will be cloned and will have new unique IDs. diff --git a/openmc/surface.py b/openmc/surface.py index 2ffbf0722..55e3f459d 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -2209,6 +2209,21 @@ class Halfspace(Region): surfaces[self.surface.id] = self.surface return surfaces + def remove_redundant_surfaces(self, redundant_surfaces): + """Recursively remove all redundant surfaces referenced by this region + + Parameters + ---------- + redundant_surfaces : dict + Dictionary mapping redundant surface IDs to surface IDs for the + :class:`openmc.Surface` instances that should replace them. + + """ + + surf = redundant_surfaces.get(self.surface.id) + if surf is not None: + self.surface = surf + def clone(self, memo=None): """Create a copy of this halfspace, with a cloned surface with a unique ID. diff --git a/src/bank.cpp b/src/bank.cpp index 172ecf84f..e91a1c1b9 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -1,12 +1,11 @@ #include "openmc/bank.h" - #include "openmc/capi.h" #include "openmc/error.h" +#include "openmc/simulation.h" +#include "openmc/message_passing.h" #include -// Explicit template instantiation definition -template class std::vector; namespace openmc { @@ -17,11 +16,21 @@ namespace openmc { namespace simulation { std::vector source_bank; -std::vector fission_bank; -std::vector secondary_bank; -#ifdef _OPENMP -std::vector master_fission_bank; -#endif + +// The fission bank is allocated as a pointer, 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 +// added to it by performing atomic incrementing captures on the fission_bank_length +// variable. A vector is not appropriate, as use of push_back() and size() +// functions would cause undefined or unintended behavior. +std::unique_ptr fission_bank; +int64_t fission_bank_length {0}; +int64_t fission_bank_max; + +// Each entry in this vector corresponds to the number of progeny produced +// this generation for the particle located at that index. This vector is +// used to efficiently sort the fission bank after each iteration. +std::vector progeny_per_particle; } // namespace simulation @@ -32,13 +41,70 @@ std::vector master_fission_bank; void free_memory_bank() { simulation::source_bank.clear(); - #pragma omp parallel - { - simulation::fission_bank.clear(); + simulation::fission_bank.reset(); + simulation::fission_bank_length = 0; + simulation::progeny_per_particle.clear(); +} + +void init_fission_bank(int64_t max) +{ + simulation::fission_bank_max = max; + simulation::fission_bank = std::make_unique(max); + simulation::fission_bank_length = 0; + simulation::progeny_per_particle.resize(simulation::work_per_rank); +} + +// Performs an O(n) sort on the fission bank, by leveraging +// the parent_id and progeny_id fields of banked particles. See the following +// paper for more details: +// "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and +// T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear +// Society, Volume 65, Page 235. +void sort_fission_bank() +{ + // Ensure we don't read off the end of the array if we ran with 0 particles + if (simulation::progeny_per_particle.size() == 0) { + return; } -#ifdef _OPENMP - simulation::master_fission_bank.clear(); -#endif + + // Perform exclusive scan summation to determine starting indices in fission + // bank for each parent particle id + int64_t tmp = simulation::progeny_per_particle[0]; + simulation::progeny_per_particle[0] = 0; + for (int64_t i = 1; i < simulation::progeny_per_particle.size(); i++) { + int64_t value = simulation::progeny_per_particle[i-1] + tmp; + tmp = simulation::progeny_per_particle[i]; + simulation::progeny_per_particle[i] = value; + } + + // TODO: C++17 introduces the exclusive_scan() function which could be + // used to replace everything above this point in this function. + + // We need a scratch vector to make permutation of the fission bank into + // sorted order easy. Under normal usage conditions, the fission bank is + // over provisioned, so we can use that as scratch space. + Particle::Bank* sorted_bank; + std::vector sorted_bank_holder; + + // If there is not enough space, allocate a temporary vector and point to it + if (simulation::fission_bank_length > simulation::fission_bank_max / 2) { + sorted_bank_holder.resize(simulation::fission_bank_length); + sorted_bank = sorted_bank_holder.data(); + } else { // otherwise, point sorted_bank to unused portion of the fission bank + sorted_bank = &simulation::fission_bank[simulation::fission_bank_length]; + } + + // Use parent and progeny indices to sort fission bank + for (int64_t i = 0; i < simulation::fission_bank_length; i++) { + const auto& site = simulation::fission_bank[i]; + int64_t offset = site.parent_id - 1 - simulation::work_index[mpi::rank]; + int64_t idx = simulation::progeny_per_particle[offset] + site.progeny_id; + sorted_bank[idx] = site; + } + + // Copy sorted bank into the fission bank + std::copy(sorted_bank, sorted_bank + simulation::fission_bank_length, + simulation::fission_bank.get()); } //============================================================================== @@ -59,12 +125,12 @@ extern "C" int openmc_source_bank(void** ptr, int64_t* n) extern "C" int openmc_fission_bank(void** ptr, int64_t* n) { - if (simulation::fission_bank.size() == 0) { + if (simulation::fission_bank_length == 0) { set_errmsg("Fission bank has not been allocated."); return OPENMC_E_ALLOCATE; } else { - *ptr = simulation::fission_bank.data(); - *n = simulation::fission_bank.size(); + *ptr = simulation::fission_bank.get(); + *n = simulation::fission_bank_length; return 0; } } diff --git a/src/cell.cpp b/src/cell.cpp index b0b074282..d56e877d3 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -500,7 +500,7 @@ CSGCell::contains(Position r, Direction u, int32_t on_surface) const //============================================================================== std::pair -CSGCell::distance(Position r, Direction u, int32_t on_surface) const +CSGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; @@ -786,13 +786,13 @@ CSGCell::contains_complex(Position r, Direction u, int32_t on_surface) const DAGCell::DAGCell() : Cell{} {}; std::pair -DAGCell::distance(Position r, Direction u, int32_t on_surface) const +DAGCell::distance(Position r, Direction u, int32_t on_surface, Particle* p) const { // if we've changed direction or we're not on a surface, // reset the history and update last direction - if (u != simulation::last_dir || on_surface == 0) { - simulation::history.reset(); - simulation::last_dir = u; + if (u != p->last_dir_ || on_surface == 0) { + p->history_.reset(); + p->last_dir_ = u; } moab::ErrorCode rval; @@ -801,7 +801,7 @@ DAGCell::distance(Position r, Direction u, int32_t on_surface) const double dist; double pnt[3] = {r.x, r.y, r.z}; double dir[3] = {u.x, u.y, u.z}; - rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &simulation::history); + rval = dagmc_ptr_->ray_fire(vol, pnt, dir, hit_surf, dist, &p->history_); MB_CHK_ERR_CONT(rval); int surf_idx; if (hit_surf != 0) { diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 33b2e8af8..e14a6062a 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -39,13 +39,6 @@ namespace openmc { const std::string DAGMC_FILENAME = "dagmc.h5m"; -namespace simulation { - -moab::DagMC::RayHistory history; -Direction last_dir; - -} - namespace model { moab::DagMC* DAG; diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 1fc1d4ca9..35d42ba1d 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -20,10 +20,6 @@ #include "openmc/timer.h" #include "openmc/tallies/tally.h" -#ifdef _OPENMP -#include -#endif - #include // for min #include #include // for sqrt, abs, pow @@ -87,7 +83,7 @@ void synchronize_bank() #ifdef OPENMC_MPI int64_t start = 0; - int64_t n_bank = simulation::fission_bank.size(); + int64_t n_bank = simulation::fission_bank_length; MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); // While we would expect the value of start on rank 0 to be 0, the MPI @@ -95,14 +91,14 @@ void synchronize_bank() // significant if (mpi::rank == 0) start = 0; - int64_t finish = start + simulation::fission_bank.size(); + int64_t finish = start + simulation::fission_bank_length; int64_t total = finish; MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm); #else int64_t start = 0; - int64_t finish = simulation::fission_bank.size(); - int64_t total = simulation::fission_bank.size(); + int64_t finish = simulation::fission_bank_length; + int64_t total = finish; #endif // If there are not that many particles per generation, it's possible that no @@ -110,7 +106,7 @@ void synchronize_bank() // extra logic to treat this circumstance, we really want to ensure the user // runs enough particles to avoid this in the first place. - if (simulation::fission_bank.empty()) { + if (simulation::fission_bank_length == 0) { fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank)); } @@ -142,7 +138,9 @@ void synchronize_bank() int64_t index_temp = 0; std::vector temp_sites(3*simulation::work_per_rank); - for (const auto& site : simulation::fission_bank) { + for (int64_t i = 0; i < simulation::fission_bank_length; i++ ) { + const auto& site = simulation::fission_bank[i]; + // If there are less than n_particles particles banked, automatically add // int(n_particles/total) sites to temp_sites. For example, if you need // 1000 and 300 were banked, this would add 3 source sites per banked site @@ -197,7 +195,7 @@ void synchronize_bank() // fission bank sites_needed = settings::n_particles - finish; for (int i = 0; i < sites_needed; ++i) { - int i_bank = simulation::fission_bank.size() - sites_needed + i; + int i_bank = simulation::fission_bank_length - sites_needed + i; temp_sites[index_temp] = simulation::fission_bank[i_bank]; ++index_temp; } @@ -351,40 +349,6 @@ void calculate_average_keff() } } -#ifdef _OPENMP -void join_bank_from_threads() -{ - int n_threads = omp_get_max_threads(); - - #pragma omp parallel - { - // Copy thread fission bank sites to one shared copy - #pragma omp for ordered schedule(static) - for (int i = 0; i < n_threads; ++i) { - #pragma omp ordered - { - std::copy( - simulation::fission_bank.cbegin(), - simulation::fission_bank.cend(), - std::back_inserter(simulation::master_fission_bank) - ); - } - } - - // Make sure all threads have made it to this point - #pragma omp barrier - - // Now copy the shared fission bank sites back to the master thread's copy. - if (omp_get_thread_num() == 0) { - simulation::fission_bank = simulation::master_fission_bank; - simulation::master_fission_bank.clear(); - } else { - simulation::fission_bank.clear(); - } - } -} -#endif - int openmc_get_keff(double* k_combined) { k_combined[0] = 0.0; @@ -542,7 +506,8 @@ void shannon_entropy() // Get source weight in each mesh bin bool sites_outside; xt::xtensor p = simulation::entropy_mesh->count_sites( - simulation::fission_bank, &sites_outside); + simulation::fission_bank.get(), simulation::fission_bank_length, + &sites_outside); // display warning message if there were sites outside entropy box if (sites_outside) { @@ -581,7 +546,7 @@ void ufs_count_sites() // count number of source sites in each ufs mesh cell bool sites_outside; simulation::source_frac = simulation::ufs_mesh->count_sites( - simulation::source_bank, &sites_outside); + simulation::source_bank.data(), simulation::source_bank.size(), &sites_outside); // Check for sites outside of the mesh if (mpi::master && sites_outside) { diff --git a/src/geometry.cpp b/src/geometry.cpp index 48a75af46..8bbc4d90c 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -119,7 +119,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } // Announce the cell that the particle is entering. - if (found && (settings::verbosity >= 10 || simulation::trace)) { + if (found && (settings::verbosity >= 10 || p->trace_)) { std::stringstream msg; msg << " Entering cell " << model::cells[i_cell]->id_; write_message(msg, 1); @@ -297,7 +297,7 @@ cross_lattice(Particle* p, const BoundaryInfo& boundary) auto& coord {p->coord_[p->n_coord_ - 1]}; auto& lat {*model::lattices[coord.lattice]}; - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || p->trace_) { std::stringstream msg; msg << " Crossing lattice " << lat.id_ << ". Current position (" << coord.lattice_x << "," << coord.lattice_y << "," @@ -370,7 +370,7 @@ BoundaryInfo distance_to_boundary(Particle* p) Cell& c {*model::cells[coord.cell]}; // Find the oncoming surface in this cell and the distance to it. - auto surface_distance = c.distance(r, u, p->surface_); + auto surface_distance = c.distance(r, u, p->surface_, p); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; diff --git a/src/initialize.cpp b/src/initialize.cpp index 38b8e5c29..166fddcf8 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -101,18 +101,20 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype Particle::Bank b; - MPI_Aint disp[6]; + MPI_Aint disp[8]; 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]); - for (int i = 5; i >= 0; --i) disp[i] -= disp[0]; + MPI_Get_address(&b.parent_id, &disp[6]); + MPI_Get_address(&b.progeny_id, &disp[7]); + for (int i = 7; i >= 0; --i) disp[i] -= disp[0]; - int blocks[] {3, 3, 1, 1, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT}; - MPI_Type_create_struct(6, blocks, disp, types, &mpi::bank); + 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); MPI_Type_commit(&mpi::bank); } #endif // OPENMC_MPI diff --git a/src/mesh.cpp b/src/mesh.cpp index f679c402a..06487bea6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -790,7 +790,7 @@ void RegularMesh::to_hdf5(hid_t group) const } xt::xtensor -RegularMesh::count_sites(const std::vector& bank, +RegularMesh::count_sites(const Particle::Bank* bank, int64_t length, bool* outside) const { // Determine shape of array for counts @@ -801,7 +801,9 @@ RegularMesh::count_sites(const std::vector& bank, xt::xarray cnt {shape, 0.0}; bool outside_ = false; - for (const auto& site : bank) { + for (int64_t i = 0; i < length; i++) { + const auto& site = bank[i]; + // determine scoring bin for entropy mesh int mesh_bin = get_bin(site.r); diff --git a/src/output.cpp b/src/output.cpp index faa9a7acd..220f3cff4 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -474,6 +474,8 @@ void print_runtime() show_time("Time synchronizing fission bank", time_bank.elapsed(), 1); show_time("Sampling source sites", time_bank_sample.elapsed(), 2); show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2); + } else { + show_time("Time sampling source", time_sample_source.elapsed(), 1); } show_time("Time accumulating tallies", time_tallies.elapsed(), 1); show_time("Total time for finalization", time_finalize.elapsed()); @@ -678,9 +680,14 @@ write_tallies() } } + // Initialize Filter Matches Object + std::vector filter_matches; + // Allocate space for tally filter matches + filter_matches.resize(model::tally_filters.size()); + // Loop over all filter bin combinations. - auto filter_iter = FilterBinIter(tally, false); - auto end = FilterBinIter(tally, true); + auto filter_iter = FilterBinIter(tally, false, &filter_matches); + auto end = FilterBinIter(tally, true, &filter_matches); for (; filter_iter != end; ++filter_iter) { auto filter_index = filter_iter.index_; @@ -691,7 +698,7 @@ write_tallies() if (filter_index % tally.strides(i) == 0) { auto i_filt = tally.filters(i); const auto& filt {*model::tally_filters[i_filt]}; - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches[i_filt]}; tallies_out << std::string(indent+1, ' ') << filt.text_label(match.i_bin_) << "\n"; } diff --git a/src/particle.cpp b/src/particle.cpp index 8bcde89d2..4e836a8dc 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -84,9 +84,9 @@ Particle::clear() void Particle::create_secondary(Direction u, double E, Type type) { - simulation::secondary_bank.emplace_back(); + secondary_bank_.emplace_back(); - auto& bank {simulation::secondary_bank.back()}; + auto& bank {secondary_bank_.back()}; bank.particle = type; bank.wgt = wgt_; bank.r = this->r(); @@ -107,6 +107,7 @@ Particle::from_source(const Bank* src) material_ = C_NONE; n_collision_ = 0; fission_ = false; + std::fill(flux_derivs_.begin(), flux_derivs_.end(), 0.0); // Copy attributes from source bank site type_ = src->particle; @@ -129,258 +130,244 @@ Particle::from_source(const Bank* src) } void -Particle::transport() +Particle::event_calculate_xs() { - // Display message if high verbosity or trace is on - if (settings::verbosity >= 9 || simulation::trace) { - write_message("Simulating Particle " + std::to_string(id_)); + // Set the random number stream + if (type_ == Particle::Type::neutron) { + stream_ = STREAM_TRACKING; + } else { + stream_ = STREAM_PHOTON; } - // Initialize number of events to zero - int n_event = 0; + // Store pre-collision particle properties + wgt_last_ = wgt_; + E_last_ = E_; + u_last_ = this->u(); + r_last_ = this->r(); - // Add paricle's starting weight to count for normalizing tallies later - #pragma omp atomic - simulation::total_weight += wgt_; + // Reset event variables + event_ = EVENT_KILL; + event_nuclide_ = NUCLIDE_NONE; + event_mt_ = REACTION_NONE; + + // If the cell hasn't been determined based on the particle's location, + // initiate a search for the current cell. This generally happens at the + // beginning of the history and again for any secondary particles + if (coord_[n_coord_ - 1].cell == C_NONE) { + if (!find_cell(this, false)) { + this->mark_as_lost("Could not find the cell containing particle " + + std::to_string(id_)); + return; + } + + // Set birth cell attribute + if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; + } + + // Write particle track. + if (write_track_) write_particle_track(*this); + + if (settings::check_overlaps) check_cell_overlap(this); + + // Calculate microscopic and macroscopic cross sections + if (material_ != MATERIAL_VOID) { + if (settings::run_CE) { + if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) { + // If the material is the same as the last material and the + // temperature hasn't changed, we don't need to lookup cross + // sections again. + model::materials[material_]->calculate_xs(*this); + } + } else { + // Get the MG data; unlike the CE case above, we have to re-calculate + // cross sections for every collision since the cross sections may + // be angle-dependent + data::mg.macro_xs_[material_].calculate_xs(*this); + + // Update the particle's group while we know we are multi-group + g_last_ = g_; + } + } else { + macro_xs_.total = 0.0; + macro_xs_.absorption = 0.0; + macro_xs_.fission = 0.0; + macro_xs_.nu_fission = 0.0; + } +} + +void +Particle::event_advance() +{ + // Find the distance to the nearest boundary + boundary_ = distance_to_boundary(this); + + // Sample a distance to collision + if (type_ == Particle::Type::electron || + type_ == Particle::Type::positron) { + collision_distance_ = 0.0; + } else if (macro_xs_.total == 0.0) { + collision_distance_ = INFINITY; + } else { + collision_distance_ = -std::log(prn(this->current_seed())) / macro_xs_.total; + } + + // Select smaller of the two distances + double distance = std::min(boundary_.distance, collision_distance_); + + // Advance particle + for (int j = 0; j < n_coord_; ++j) { + coord_[j].r += distance * coord_[j].u; + } + + // Score track-length tallies + if (!model::active_tracklength_tallies.empty()) { + score_tracklength_tally(this, distance); + } + + // Score track-length estimate of k-eff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type_ == Particle::Type::neutron) { + keff_tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission; + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) { + score_track_derivative(this, distance); + } +} + +void +Particle::event_cross_surface() +{ + // Set surface that particle is on and adjust coordinate levels + surface_ = boundary_.surface_index; + n_coord_ = boundary_.coord_level; + + // Saving previous cell data + for (int j = 0; j < n_coord_; ++j) { + cell_last_[j] = coord_[j].cell; + } + n_coord_last_ = n_coord_; + + if (boundary_.lattice_translation[0] != 0 || + boundary_.lattice_translation[1] != 0 || + boundary_.lattice_translation[2] != 0) { + // Particle crosses lattice boundary + cross_lattice(this, boundary_); + event_ = EVENT_LATTICE; + } else { + // Particle crosses surface + this->cross_surface(); + event_ = EVENT_SURFACE; + } + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } +} + +void +Particle::event_collide() +{ + // Score collision estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type_ == Particle::Type::neutron) { + keff_tally_collision_ += wgt_ * macro_xs_.nu_fission + / macro_xs_.total; + } + + // Score surface current tallies -- this has to be done before the collision + // since the direction of the particle will change and we need to use the + // pre-collision direction to figure out what mesh surfaces were crossed + + if (!model::active_meshsurf_tallies.empty()) + score_surface_tally(this, model::active_meshsurf_tallies); + + // Clear surface component + surface_ = 0; - // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { - for (auto& micro : neutron_xs_) micro.last_E = 0.0; + collision(this); + } else { + collision_mg(this); } - // Prepare to write out particle track. - if (write_track_) add_particle_track(); - - // Every particle starts with no accumulated flux derivative. - if (!model::active_tallies.empty()) zero_flux_derivs(); - - while (true) { - // Set the random number stream - if (type_ == Particle::Type::neutron) { - stream_ = STREAM_TRACKING; + // Score collision estimator tallies -- this is done after a collision + // has occurred rather than before because we need information on the + // outgoing energy for any tallies with an outgoing energy filter + if (!model::active_collision_tallies.empty()) score_collision_tally(this); + if (!model::active_analog_tallies.empty()) { + if (settings::run_CE) { + score_analog_tally_ce(this); } else { - stream_ = STREAM_PHOTON; - } - - // Store pre-collision particle properties - wgt_last_ = wgt_; - E_last_ = E_; - u_last_ = this->u(); - r_last_ = this->r(); - - // Reset event variables - event_ = TallyEvent::KILL; - event_nuclide_ = NUCLIDE_NONE; - event_mt_ = REACTION_NONE; - - // If the cell hasn't been determined based on the particle's location, - // initiate a search for the current cell. This generally happens at the - // beginning of the history and again for any secondary particles - if (coord_[n_coord_ - 1].cell == C_NONE) { - if (!find_cell(this, false)) { - this->mark_as_lost("Could not find the cell containing particle " - + std::to_string(id_)); - return; - } - - // Set birth cell attribute - if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; - } - - // Write particle track. - if (write_track_) write_particle_track(*this); - - if (settings::check_overlaps) check_cell_overlap(this); - - // Calculate microscopic and macroscopic cross sections - if (material_ != MATERIAL_VOID) { - if (settings::run_CE) { - if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) { - // If the material is the same as the last material and the - // temperature hasn't changed, we don't need to lookup cross - // sections again. - model::materials[material_]->calculate_xs(*this); - } - } else { - // Get the MG data; unlike the CE case above, we have to re-calculate - // cross sections for every collision since the cross sections may - // be angle-dependent - data::mg.macro_xs_[material_].calculate_xs(*this); - - // Update the particle's group while we know we are multi-group - g_last_ = g_; - } - } else { - macro_xs_.total = 0.0; - macro_xs_.absorption = 0.0; - macro_xs_.fission = 0.0; - macro_xs_.nu_fission = 0.0; - } - - // Find the distance to the nearest boundary - auto boundary = distance_to_boundary(this); - - // Sample a distance to collision - double d_collision; - if (type_ == Particle::Type::electron || - type_ == Particle::Type::positron) { - d_collision = 0.0; - } else if (macro_xs_.total == 0.0) { - d_collision = INFINITY; - } else { - d_collision = -std::log(prn(this->current_seed())) / macro_xs_.total; - } - - // Select smaller of the two distances - double distance = std::min(boundary.distance, d_collision); - - // Advance particle - for (int j = 0; j < n_coord_; ++j) { - coord_[j].r += distance * coord_[j].u; - } - - // Score track-length tallies - if (!model::active_tracklength_tallies.empty()) { - score_tracklength_tally(this, distance); - } - - // Score track-length estimate of k-eff - if (settings::run_mode == RunMode::EIGENVALUE && - type_ == Particle::Type::neutron) { - global_tally_tracklength += wgt_ * distance * macro_xs_.nu_fission; - } - - // Score flux derivative accumulators for differential tallies. - if (!model::active_tallies.empty()) { - score_track_derivative(this, distance); - } - - if (d_collision > boundary.distance) { - // ==================================================================== - // PARTICLE CROSSES SURFACE - - // Set surface that particle is on and adjust coordinate levels - surface_ = boundary.surface_index; - n_coord_ = boundary.coord_level; - - // Saving previous cell data - for (int j = 0; j < n_coord_; ++j) { - cell_last_[j] = coord_[j].cell; - } - n_coord_last_ = n_coord_; - - if (boundary.lattice_translation[0] != 0 || - boundary.lattice_translation[1] != 0 || - boundary.lattice_translation[2] != 0) { - // Particle crosses lattice boundary - cross_lattice(this, boundary); - event_ = TallyEvent::LATTICE; - } else { - // Particle crosses surface - this->cross_surface(); - event_ = TallyEvent::SURFACE; - } - // Score cell to cell partial currents - if (!model::active_surface_tallies.empty()) { - score_surface_tally(this, model::active_surface_tallies); - } - } else { - // ==================================================================== - // PARTICLE HAS COLLISION - - // Score collision estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE && - type_ == Particle::Type::neutron) { - global_tally_collision += wgt_ * macro_xs_.nu_fission - / macro_xs_.total; - } - - // Score surface current tallies -- this has to be done before the collision - // since the direction of the particle will change and we need to use the - // pre-collision direction to figure out what mesh surfaces were crossed - - if (!model::active_meshsurf_tallies.empty()) - score_surface_tally(this, model::active_meshsurf_tallies); - - // Clear surface component - surface_ = 0; - - if (settings::run_CE) { - collision(this); - } else { - collision_mg(this); - } - - // Score collision estimator tallies -- this is done after a collision - // has occurred rather than before because we need information on the - // outgoing energy for any tallies with an outgoing energy filter - if (!model::active_collision_tallies.empty()) score_collision_tally(this); - if (!model::active_analog_tallies.empty()) { - if (settings::run_CE) { - score_analog_tally_ce(this); - } else { - score_analog_tally_mg(this); - } - } - - // Reset banked weight during collision - n_bank_ = 0; - n_bank_second_ = 0; - wgt_bank_ = 0.0; - for (int& v : n_delayed_bank_) v = 0; - - // Reset fission logical - fission_ = false; - - // Save coordinates for tallying purposes - r_last_current_ = this->r(); - - // Set last material to none since cross sections will need to be - // re-evaluated - material_last_ = C_NONE; - - // Set all directions to base level -- right now, after a collision, only - // the base level directions are changed - for (int j = 0; j < n_coord_ - 1; ++j) { - if (coord_[j + 1].rotated) { - // If next level is rotated, apply rotation matrix - const auto& m {model::cells[coord_[j].cell]->rotation_}; - const auto& u {coord_[j].u}; - coord_[j + 1].u = u.rotate(m); - } else { - // Otherwise, copy this level's direction - coord_[j+1].u = coord_[j].u; - } - } - - // Score flux derivative accumulators for differential tallies. - if (!model::active_tallies.empty()) score_collision_derivative(this); - } - - // If particle has too many events, display warning and kill it - ++n_event; - if (n_event == MAX_EVENTS) { - warning("Particle " + std::to_string(id_) + - " underwent maximum number of events."); - alive_ = false; - } - - // Check for secondary particles if this particle is dead - if (!alive_) { - // If no secondary particles, break out of event loop - if (simulation::secondary_bank.empty()) break; - - this->from_source(&simulation::secondary_bank.back()); - simulation::secondary_bank.pop_back(); - n_event = 0; - - // Enter new particle in particle track file - if (write_track_) add_particle_track(); + score_analog_tally_mg(this); } } + // Reset banked weight during collision + n_bank_ = 0; + n_bank_second_ = 0; + wgt_bank_ = 0.0; + for (int& v : n_delayed_bank_) v = 0; + + // Reset fission logical + fission_ = false; + + // Save coordinates for tallying purposes + r_last_current_ = this->r(); + + // Set last material to none since cross sections will need to be + // re-evaluated + material_last_ = C_NONE; + + // Set all directions to base level -- right now, after a collision, only + // the base level directions are changed + for (int j = 0; j < n_coord_ - 1; ++j) { + if (coord_[j + 1].rotated) { + // If next level is rotated, apply rotation matrix + const auto& m {model::cells[coord_[j].cell]->rotation_}; + const auto& u {coord_[j].u}; + coord_[j + 1].u = u.rotate(m); + } else { + // Otherwise, copy this level's direction + coord_[j+1].u = coord_[j].u; + } + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) score_collision_derivative(this); +} + +void +Particle::event_revive_from_secondary() +{ + // If particle has too many events, display warning and kill it + ++n_event_; + if (n_event_ == MAX_EVENTS) { + warning("Particle " + std::to_string(id_) + + " underwent maximum number of events."); + alive_ = false; + } + + // Check for secondary particles if this particle is dead + if (!alive_) { + // If no secondary particles, break out of event loop + if (secondary_bank_.empty()) return; + + this->from_source(&secondary_bank_.back()); + secondary_bank_.pop_back(); + n_event_ = 0; + + // Enter new particle in particle track file + if (write_track_) add_particle_track(*this); + } +} + +void +Particle::event_death() +{ #ifdef DAGMC - if (settings::dagmc) simulation::history.reset(); + if (settings::dagmc) history_.reset(); #endif // Finish particle track output. @@ -388,15 +375,39 @@ Particle::transport() write_particle_track(*this); finalize_particle_track(*this); } + + // Contribute tally reduction variables to global accumulator + #pragma omp atomic + global_tally_absorption += keff_tally_absorption_; + #pragma omp atomic + global_tally_collision += keff_tally_collision_; + #pragma omp atomic + global_tally_tracklength += keff_tally_tracklength_; + #pragma omp atomic + global_tally_leakage += keff_tally_leakage_; + + // Reset particle tallies once accumulated + keff_tally_absorption_ = 0.0; + keff_tally_collision_ = 0.0; + keff_tally_tracklength_ = 0.0; + keff_tally_leakage_ = 0.0; + + // Record the number of progeny created by this particle. + // This data will be used to efficiently sort the fission bank. + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + int64_t offset = id_ - 1 - simulation::work_index[mpi::rank]; + simulation::progeny_per_particle[offset] = n_progeny_; + } } + void Particle::cross_surface() { int i_surface = std::abs(surface_); // TODO: off-by-one const auto& surf {model::surfaces[i_surface - 1].get()}; - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || trace_) { write_message(" Crossing surface " + std::to_string(surf->id_)); } @@ -420,10 +431,10 @@ Particle::cross_surface() } // Score to global leakage tally - global_tally_leakage += wgt_; + keff_tally_leakage_ += wgt_; // Display message - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || trace_) { write_message(" Leaked out of surface " + std::to_string(surf->id_)); } return; @@ -460,7 +471,7 @@ Particle::cross_surface() } Direction u = (surf->bc_ == Surface::BoundaryType::REFLECT) ? - surf->reflect(this->r(), this->u()) : + surf->reflect(this->r(), this->u(), this) : surf->diffuse_reflect(this->r(), this->u(), this->current_seed()); // Make sure new particle direction is normalized @@ -487,7 +498,7 @@ Particle::cross_surface() r_last_current_ = this->r() + TINY_BIT*this->u(); // Diagnostic message - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || trace_) { write_message(" Reflected from surface " + std::to_string(surf->id_)); } return; @@ -541,7 +552,7 @@ Particle::cross_surface() r_last_current_ = this->r() + TINY_BIT*this->u(); // Diagnostic message - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || trace_) { write_message(" Hit periodic boundary on surface " + std::to_string(surf->id_)); } @@ -667,7 +678,7 @@ Particle::write_restart() const write_dataset(file_id, "id", id_); write_dataset(file_id, "type", static_cast(type_)); - int64_t i = simulation::current_work; + int64_t i = current_work_; write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt); write_dataset(file_id, "energy", simulation::source_bank[i-1].E); write_dataset(file_id, "xyz", simulation::source_bank[i-1].r); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 1f5ff9eb8..97b82c5c6 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -10,7 +10,9 @@ #include "openmc/random_lcg.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/tallies/derivative.h" #include "openmc/tallies/tally.h" +#include "openmc/track_output.h" #include // for copy #include @@ -100,8 +102,25 @@ void run_particle_restart() } init_particle_seeds(particle_seed, p.seeds_); + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (auto& micro : p.neutron_xs_) micro.last_E = 0.0; + } + + // Prepare to write out particle track. + if (p.write_track_) add_particle_track(p); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) { + p.flux_derivs_.resize(model::tally_derivs.size(), 0.0); + std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0); + } + + // Allocate space for tally filter matches + p.filter_matches_.resize(model::tally_filters.size()); + // Transport neutron - p.transport(); + transport_history_based_single_particle(p); // Write output if particle made it print_particle(&p); diff --git a/src/physics.cpp b/src/physics.cpp index 483424493..b6a509bb0 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -60,7 +60,7 @@ void collision(Particle* p) } // Display information about collision - if (settings::verbosity >= 10 || simulation::trace) { + if (settings::verbosity >= 10 || p->trace_) { std::stringstream msg; if (p->event_ == TallyEvent::KILL) { msg << " Killed. Energy = " << p->E_ << " eV."; @@ -95,14 +95,14 @@ void sample_neutron_reaction(Particle* p) if (nuc->fissionable_) { Reaction* rx = sample_fission(i_nuclide, p); if (settings::run_mode == RunMode::EIGENVALUE) { - create_fission_sites(p, i_nuclide, rx, simulation::fission_bank); + create_fission_sites(p, i_nuclide, rx); } else if (settings::run_mode == RunMode::FIXED_SOURCE && settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank); + create_fission_sites(p, i_nuclide, rx); // Make sure particle population doesn't grow out of control for // subcritical multiplication problems. - if (simulation::secondary_bank.size() >= 10000) { + if (p->secondary_bank_.size() >= 10000) { fatal_error("The secondary particle bank appears to be growing without " "bound. You are likely running a subcritical multiplication problem " "with k-effective close to or greater than one."); @@ -146,8 +146,7 @@ void sample_neutron_reaction(Particle* p) } void -create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, - std::vector& bank) +create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx) { // If uniform fission source weighting is turned on, we increase or decrease // the expected number of fission sites produced @@ -169,28 +168,74 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; - p->fission_ = true; - for (int i = 0; i < nu; ++i) { - // Create new bank site and get reference to last element - bank.emplace_back(); - auto& site {bank.back()}; + // Clear out particle's nu fission bank + p->nu_bank_.clear(); - // Bank source neutrons by copying the particle data - site.r = p->r(); - site.particle = Particle::Type::neutron; - site.wgt = 1. / weight; + p->fission_ = true; + int skipped = 0; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE); + + for (int i = 0; i < nu; ++i) { + Particle::Bank* site; + if (use_fission_bank) { + int64_t idx; + #pragma omp atomic capture + idx = simulation::fission_bank_length++; + if (idx >= simulation::fission_bank_max) { + warning("The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked."); + #pragma omp atomic write + simulation::fission_bank_length = simulation::fission_bank_max; + skipped++; + break; + } + site = &simulation::fission_bank[idx]; + } else { + // Create new bank site and get reference to last element + auto& bank = p->secondary_bank_; + bank.emplace_back(); + site = &bank.back(); + } + site->r = p->r(); + site->particle = Particle::Type::neutron; + site->wgt = 1. / weight; + site->parent_id = p->id_; + site->progeny_id = p->n_progeny_++; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->current_seed()); + sample_fission_neutron(i_nuclide, rx, p->E_, site, p->current_seed()); // Set the delayed group on the particle as well - p->delayed_group_ = site.delayed_group; + p->delayed_group_ = site->delayed_group; // Increment the number of neutrons born delayed if (p->delayed_group_ > 0) { nu_d[p->delayed_group_-1]++; } + + // Write fission particles to nuBank + if (use_fission_bank) { + p->nu_bank_.emplace_back(); + Particle::NuBank* nu_bank_entry = &p->nu_bank_.back(); + nu_bank_entry->wgt = site->wgt; + nu_bank_entry->E = site->E; + nu_bank_entry->delayed_group = site->delayed_group; + } } + + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (nu == skipped) { + p->fission_ = false; + return; + } + + // If shared fission bank was full, but some fissions could be added, + // reduce nu accordingly + nu -= skipped; // Store the total weight banked for analog fission tallies p->n_bank_ = nu; @@ -548,7 +593,7 @@ void absorption(Particle* p, int i_nuclide) // Score implicit absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - global_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[ + p->keff_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[ i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption; } } else { @@ -557,7 +602,7 @@ void absorption(Particle* p, int i_nuclide) prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - global_tally_absorption += p->wgt_ * p->neutron_xs_[ + p->keff_tally_absorption_ += p->wgt_ * p->neutron_xs_[ i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption; } @@ -696,12 +741,12 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Find speed of neutron in CM vel = v_n.norm(); - // Sample scattering angle, checking if it is an ncorrelated angle-energy - // distribution + // Sample scattering angle, checking if angle distribution is present (assume + // isotropic otherwise) double mu_cm; auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); - if (d_) { + if (!d_->angle().empty()) { mu_cm = d_->angle().sample(p->E_, p->current_seed()); } else { mu_cm = 2.0*prn(p->current_seed()) - 1.0; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index efc500bee..54ab4ec4b 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -31,7 +31,7 @@ collision_mg(Particle* p) sample_reaction(p); // Display information about collision - if ((settings::verbosity >= 10) || (simulation::trace)) { + if ((settings::verbosity >= 10) || (p->trace_)) { std::stringstream msg; msg << " Energy Group = " << p->g_; write_message(msg, 1); @@ -47,11 +47,10 @@ sample_reaction(Particle* p) // absorption (including fission) if (model::materials[p->material_]->fissionable_) { - if (settings::run_mode == RunMode::EIGENVALUE) { - create_fission_sites(p, simulation::fission_bank); - } else if ((settings::run_mode == RunMode::FIXED_SOURCE) && - (settings::create_fission_neutrons)) { - create_fission_sites(p, simulation::secondary_bank); + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::create_fission_neutrons)) { + create_fission_sites(p); } } @@ -91,7 +90,7 @@ scatter(Particle* p) } void -create_fission_sites(Particle* p, std::vector& bank) +create_fission_sites(Particle* p) { // If uniform fission source weighting is turned on, we increase or decrease // the expected number of fission sites produced @@ -114,17 +113,45 @@ create_fission_sites(Particle* p, std::vector& bank) // Initialize the counter of delayed neutrons encountered for each delayed // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + // Clear out particle's nu fission bank + p->nu_bank_.clear(); p->fission_ = true; + int skipped = 0; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RUN_MODE_EIGENVALUE); + for (int i = 0; i < nu; ++i) { - // Create new bank site and get reference to last element - bank.emplace_back(); - auto& site {bank.back()}; + Particle::Bank* site; + if (use_fission_bank) { + int64_t idx; + #pragma omp atomic capture + idx = simulation::fission_bank_length++; + if (idx >= simulation::fission_bank_max) { + warning("The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked."); + #pragma omp atomic write + simulation::fission_bank_length = simulation::fission_bank_max; + skipped++; + break; + } + site = &simulation::fission_bank[idx]; + } else { + // Create new bank site and get reference to last element + auto& bank = p->secondary_bank_; + bank.emplace_back(); + site = &bank.back(); + } // Bank source neutrons by copying the particle data - site.r = p->r(); - site.particle = Particle::Type::neutron; - site.wgt = 1. / weight; + site->r = p->r(); + site->particle = Particle::Type::neutron; + site->wgt = 1. / weight; + site->parent_id = p->id_; + site->progeny_id = p->n_progeny_++; // Sample the cosine of the angle, assuming fission neutrons are emitted // isotropically @@ -132,9 +159,9 @@ create_fission_sites(Particle* p, std::vector& bank) // Sample the azimuthal angle uniformly in [0, 2.pi) double phi = 2. * PI * prn(p->current_seed() ); - site.u.x = mu; - site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); - site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); + site->u.x = mu; + site->u.y = std::sqrt(1. - mu * mu) * std::cos(phi); + site->u.z = std::sqrt(1. - mu * mu) * std::sin(phi); // Sample secondary energy distribution for the fission reaction int dg; @@ -142,10 +169,10 @@ create_fission_sites(Particle* p, std::vector& bank) data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout, p->current_seed()); // Store the energy and delayed groups on the fission bank - site.E = gout; + site->E = gout; // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest // of the code, 0 is prompt. - site.delayed_group = dg + 1; + site->delayed_group = dg + 1; // Set the delayed group on the particle as well p->delayed_group_ = dg + 1; @@ -154,7 +181,27 @@ create_fission_sites(Particle* p, std::vector& bank) if (p->delayed_group_ > 0) { nu_d[dg]++; } + + // Write fission particles to nuBank + if (use_fission_bank) { + p->nu_bank_.emplace_back(); + Particle::NuBank* nu_bank_entry = &p->nu_bank_.back(); + nu_bank_entry->wgt = site->wgt; + nu_bank_entry->E = site->E; + nu_bank_entry->delayed_group = site->delayed_group; + } } + + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (nu == skipped) { + p->fission_ = false; + return; + } + + // If shared fission bank was full, but some fissions could be added, + // reduce nu accordingly + nu -= skipped; // Store the total weight banked for analog fission tallies p->n_bank_ = nu; @@ -176,13 +223,11 @@ absorption(Particle* p) p->wgt_last_ = p->wgt_; // Score implicit absorpion estimate of keff - #pragma omp atomic - global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission / + p->keff_tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission / p->macro_xs_.absorption; } else { if (p->macro_xs_.absorption > prn(p->current_seed()) * p->macro_xs_.total) { - #pragma omp atomic - global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission / + p->keff_tally_absorption_ += p->wgt_ * p->macro_xs_.nu_fission / p->macro_xs_.absorption; p->alive_ = false; p->event_ = TallyEvent::ABSORB; diff --git a/src/simulation.cpp b/src/simulation.cpp index e2e730dc6..b3885b456 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -17,18 +17,25 @@ #include "openmc/source.h" #include "openmc/state_point.h" #include "openmc/timer.h" +#include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" #include "openmc/tallies/trigger.h" +#include "openmc/track_output.h" #ifdef _OPENMP #include #endif #include "xtensor/xview.hpp" +#ifdef OPENMC_MPI +#include +#endif + #include #include + //============================================================================== // C API functions //============================================================================== @@ -61,12 +68,6 @@ int openmc_simulation_init() // Determine how much work each process should do calculate_work(); - // Allocate array for matching filter bins - #pragma omp parallel - { - simulation::filter_matches.resize(model::tally_filters.size()); - } - // Allocate source bank, and for eigenvalue simulations also allocate the // fission bank allocate_banks(); @@ -139,11 +140,6 @@ int openmc_simulation_finalize() // Write tally results to tallies.out if (settings::output_tallies && mpi::master) write_tallies(); - #pragma omp parallel - { - simulation::filter_matches.clear(); - } - // Deactivate all tallies for (auto& t : model::tallies) { t->active_ = false; @@ -186,20 +182,8 @@ int openmc_next_batch(int* status) // Start timer for transport simulation::time_transport.start(); - // ==================================================================== - // LOOP OVER PARTICLES - - #pragma omp parallel for schedule(runtime) - for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { - simulation::current_work = i_work; - - // grab source particle from bank - Particle p; - initialize_history(&p, simulation::current_work); - - // transport particle - p.transport(); - } + // Transport loop + transport_history_based(); // Accumulate time for transport simulation::time_transport.stop(); @@ -242,7 +226,6 @@ namespace simulation { int current_batch; int current_gen; -int64_t current_work; bool initialized {false}; double keff {1.0}; double keff_std; @@ -264,8 +247,6 @@ const RegularMesh* ufs_mesh {nullptr}; std::vector k_generation; std::vector work_index; -// Threadprivate variables -bool trace; //!< flag to show debug information } // namespace simulation @@ -279,26 +260,9 @@ void allocate_banks() simulation::source_bank.resize(simulation::work_per_rank); if (settings::run_mode == RunMode::EIGENVALUE) { -#ifdef _OPENMP - // If OpenMP is being used, each thread needs its own private fission - // bank. Since the private fission banks need to be combined at the end of - // a generation, there is also a 'master_fission_bank' that is used to - // collect the sites from each thread. - - #pragma omp parallel - { - if (omp_get_thread_num() == 0) { - simulation::fission_bank.reserve(3*simulation::work_per_rank); - } else { - int n_threads = omp_get_num_threads(); - simulation::fission_bank.reserve(3*simulation::work_per_rank / n_threads); - } - } - simulation::master_fission_bank.reserve(3*simulation::work_per_rank); -#else - simulation::fission_bank.reserve(3*simulation::work_per_rank); -#endif + init_fission_bank(3*simulation::work_per_rank); } + } void initialize_batch() @@ -399,7 +363,7 @@ void initialize_generation() { if (settings::run_mode == RunMode::EIGENVALUE) { // Clear out the fission bank - simulation::fission_bank.clear(); + simulation::fission_bank_length = 0; // Count source sites if using uniform fission source weighting if (settings::ufs_on) ufs_count_sites(); @@ -414,34 +378,27 @@ void finalize_generation() { auto& gt = simulation::global_tallies; - // Update global tallies with the omp private accumulation variables - #pragma omp parallel - { - #pragma omp critical(increment_global_tallies) - { - if (settings::run_mode == RunMode::EIGENVALUE) { - gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision; - gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption; - gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength; - } - gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; - } - - // reset threadprivate tallies - if (settings::run_mode == RunMode::EIGENVALUE) { - global_tally_collision = 0.0; - global_tally_absorption = 0.0; - global_tally_tracklength = 0.0; - } - global_tally_leakage = 0.0; - } - - + // Update global tallies with the accumulation variables if (settings::run_mode == RunMode::EIGENVALUE) { -#ifdef _OPENMP - // Join the fission bank from each thread into one global fission bank - join_bank_from_threads(); -#endif + gt(K_COLLISION, RESULT_VALUE) += global_tally_collision; + gt(K_ABSORPTION, RESULT_VALUE) += global_tally_absorption; + gt(K_TRACKLENGTH, RESULT_VALUE) += global_tally_tracklength; + } + gt(LEAKAGE, RESULT_VALUE) += global_tally_leakage; + + // reset tallies + if (settings::run_mode == RunMode::EIGENVALUE) { + global_tally_collision = 0.0; + global_tally_absorption = 0.0; + global_tally_tracklength = 0.0; + } + global_tally_leakage = 0.0; + + if (settings::run_mode == RunMode::EIGENVALUE) { + // If using shared memory, stable sort the fission bank (by parent IDs) + // so as to allow for reproducibility regardless of which order particles + // are run in. + sort_fission_bank(); // Distribute fission bank across processors evenly synchronize_bank(); @@ -462,7 +419,9 @@ void finalize_generation() } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // For fixed-source mode, we need to sample the external source + simulation::time_sample_source.start(); fill_source_bank_fixedsource(); + simulation::time_sample_source.stop(); } } @@ -470,20 +429,24 @@ void initialize_history(Particle* p, int64_t index_source) { // set defaults p->from_source(&simulation::source_bank[index_source - 1]); + p->current_work_ = index_source; // set identifier for particle p->id_ = simulation::work_index[mpi::rank] + index_source; + // set progeny count to zero + p->n_progeny_ = 0; + // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + p->id_; init_particle_seeds(particle_seed, p->seeds_); // set particle trace - simulation::trace = false; + p->trace_ = false; if (simulation::current_batch == settings::trace_batch && simulation::current_gen == settings::trace_gen && - p->id_ == settings::trace_particle) simulation::trace = true; + p->id_ == settings::trace_particle) p->trace_ = true; // Set particle track. p->write_track_ = false; @@ -499,6 +462,38 @@ void initialize_history(Particle* p, int64_t index_source) } } } + + // Display message if high verbosity or trace is on + if (settings::verbosity >= 9 || p->trace_) { + write_message("Simulating Particle " + std::to_string(p->id_)); + } + + // Add paricle's starting weight to count for normalizing tallies later + #pragma omp atomic + simulation::total_weight += p->wgt_; + + initialize_history_partial(p); +} + +void initialize_history_partial(Particle* p) +{ + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (auto& micro : p->neutron_xs_) micro.last_E = 0.0; + } + + // Prepare to write out particle track. + if (p->write_track_) add_particle_track(*p); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) + { + p->flux_derivs_.resize(model::tally_derivs.size(), 0.0); + std::fill(p->flux_derivs_.begin(), p->flux_derivs_.end(), 0.0); + } + + // Allocate space for tally filter matches + p->filter_matches_.resize(model::tally_filters.size()); } int overall_generation() @@ -571,4 +566,31 @@ void free_memory_simulation() simulation::entropy.clear(); } +void transport_history_based_single_particle(Particle& p) +{ + while (true) { + p.event_calculate_xs(); + p.event_advance(); + if (p.collision_distance_ > p.boundary_.distance) { + p.event_cross_surface(); + } else { + p.event_collide(); + } + p.event_revive_from_secondary(); + if (!p.alive_) + break; + } + p.event_death(); +} + +void transport_history_based() +{ + #pragma omp parallel for schedule(runtime) + for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { + Particle p; + initialize_history(&p, i_work); + transport_history_based_single_particle(p); + } +} + } // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index d2e13a173..6314eb23f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -328,6 +328,7 @@ void free_memory_source() void fill_source_bank_fixedsource() { if (settings::path_source.empty()) { + #pragma omp parallel for for (int64_t i = 0; i < simulation::work_per_rank; ++i) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation()) * diff --git a/src/state_point.cpp b/src/state_point.cpp index cf733dd77..39f7a198c 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -283,6 +283,8 @@ openmc_statepoint_write(const char* filename, bool* write_source) write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed()); write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed()); write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed()); + } else { + write_dataset(runtime_group, "sampling source sites", time_sample_source.elapsed()); } write_dataset(runtime_group, "accumulating tallies", time_tallies.elapsed()); write_dataset(runtime_group, "total", time_total.elapsed()); diff --git a/src/surface.cpp b/src/surface.cpp index c1f4885be..f5399fa6a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -175,7 +175,7 @@ Surface::sense(Position r, Direction u) const } Direction -Surface::reflect(Position r, Direction u) const +Surface::reflect(Position r, Direction u, Particle* p) const { // Determine projection of direction onto normal and squared magnitude of // normal. @@ -283,11 +283,11 @@ Direction DAGSurface::normal(Position r) const return dir; } -Direction DAGSurface::reflect(Position r, Direction u) const +Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const { - simulation::history.reset_to_last_intersection(); - simulation::last_dir = Surface::reflect(r, u); - return simulation::last_dir; + p->history_.reset_to_last_intersection(); + p->last_dir_ = Surface::reflect(r, u, p); + return p->last_dir_; } void DAGSurface::to_hdf5(hid_t group_id) const {} diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 74bf122d6..1d2d2fc99 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -81,13 +81,9 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) void read_tally_derivatives(pugi::xml_node node) { - // Populate the derivatives array. This must be done in parallel because - // the derivatives are threadprivate. - #pragma omp parallel - { - for (auto deriv_node : node.children("derivative")) - model::tally_derivs.emplace_back(deriv_node); - } + // Populate the derivatives array. + for (auto deriv_node : node.children("derivative")) + model::tally_derivs.emplace_back(deriv_node); // Fill the derivative map. for (auto i = 0; i < model::tally_derivs.size(); ++i) { @@ -120,7 +116,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, // perturbated variable. const auto& deriv {model::tally_derivs[tally.deriv_]}; - auto flux_deriv = deriv.flux_deriv; + const auto flux_deriv = p->flux_derivs_[tally.deriv_]; // Handle special cases where we know that d_c/d_p must be zero. if (score_bin == SCORE_FLUX) { @@ -567,13 +563,15 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, } void -score_track_derivative(const Particle* p, double distance) +score_track_derivative(Particle* p, double distance) { // A void material cannot be perturbed so it will not affect flux derivatives. if (p->material_ == MATERIAL_VOID) return; const Material& material {*model::materials[p->material_]}; - - for (auto& deriv : model::tally_derivs) { + + for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { + const auto& deriv = model::tally_derivs[idx]; + auto& flux_deriv = p->flux_derivs_[idx]; if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -582,7 +580,7 @@ score_track_derivative(const Particle* p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist - deriv.flux_deriv -= distance * p->macro_xs_.total + flux_deriv -= distance * p->macro_xs_.total / material.density_gpcc_; break; @@ -590,7 +588,7 @@ score_track_derivative(const Particle* p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist - deriv.flux_deriv -= distance + flux_deriv -= distance * p->neutron_xs_[deriv.diff_nuclide].total; break; @@ -604,7 +602,7 @@ score_track_derivative(const Particle* p, double distance) double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_); - deriv.flux_deriv -= distance * (dsig_s + dsig_a) + flux_deriv -= distance * (dsig_s + dsig_a) * material.atom_density_(i); } } @@ -613,14 +611,17 @@ score_track_derivative(const Particle* p, double distance) } } -void score_collision_derivative(const Particle* p) +void score_collision_derivative(Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. if (p->material_ == MATERIAL_VOID) return; const Material& material {*model::materials[p->material_]}; - for (auto& deriv : model::tally_derivs) { + for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { + const auto& deriv = model::tally_derivs[idx]; + auto& flux_deriv = p->flux_derivs_[idx]; + if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -629,7 +630,7 @@ void score_collision_derivative(const Particle* p) // phi is proportional to Sigma_s // (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s // (1 / phi) * (d_phi / d_rho) = 1 / rho - deriv.flux_deriv += 1. / material.density_gpcc_; + flux_deriv += 1. / material.density_gpcc_; break; case DerivativeVariable::NUCLIDE_DENSITY: @@ -650,7 +651,7 @@ void score_collision_derivative(const Particle* p) // (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s // (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s // (1 / phi) * (d_phi / d_N) = 1 / N - deriv.flux_deriv += 1. / material.atom_density_(i); + flux_deriv += 1. / material.atom_density_(i); break; case DerivativeVariable::TEMPERATURE: @@ -665,7 +666,7 @@ void score_collision_derivative(const Particle* p) double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); - deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption); + flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption); // Note that this is an approximation! The real scattering cross // section is // Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u). @@ -680,9 +681,4 @@ void score_collision_derivative(const Particle* p) } } -void zero_flux_derivs() -{ - for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.; -} - }// namespace openmc diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index b1b6f4d61..6e416dd8d 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -39,10 +39,6 @@ namespace openmc { // Global variables //============================================================================== -namespace simulation { - std::vector filter_matches; -} - namespace model { std::vector> tally_filters; std::unordered_map filter_map; diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 20fffb737..3d5fdeff7 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1044,10 +1044,7 @@ setup_active_tallies() void free_memory_tally() { - #pragma omp parallel - { - model::tally_derivs.clear(); - } + model::tally_derivs.clear(); model::tally_deriv_map.clear(); model::tally_filters.clear(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index af4cbccaf..8ccd80261 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -26,13 +26,13 @@ namespace openmc { // FilterBinIter implementation //============================================================================== -FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p) - : tally_{tally} +FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) + : tally_{tally}, filter_matches_{p->filter_matches_} { // Find all valid bins in each relevant filter if they have not already been // found for this event. for (auto i_filt : tally_.filters()) { - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); @@ -55,8 +55,9 @@ FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p) this->compute_index_weight(); } -FilterBinIter::FilterBinIter(const Tally& tally, bool end) - : tally_{tally} +FilterBinIter::FilterBinIter(const Tally& tally, bool end, + std::vector* particle_filter_matches) + : tally_{tally}, filter_matches_{*particle_filter_matches} { // Handle the special case for an iterator that points to the end. if (end) { @@ -65,7 +66,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end) } for (auto i_filt : tally_.filters()) { - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); @@ -97,7 +98,7 @@ FilterBinIter::operator++() bool done_looping = true; for (int i = tally_.filters().size()-1; i >= 0; --i) { auto i_filt = tally_.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (match.i_bin_ < match.bins_.size()-1) { // The bin for this filter can be incremented. Increment it and do not // touch any of the remaining filters. @@ -130,7 +131,7 @@ FilterBinIter::compute_index_weight() weight_ = 1.; for (auto i = 0; i < tally_.filters().size(); ++i) { auto i_filt = tally_.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; auto i_bin = match.i_bin_; index_ += match.bins_[i_bin] * tally_.strides(i); weight_ *= match.weights_[i_bin]; @@ -144,12 +145,13 @@ FilterBinIter::compute_index_weight() //! Helper function used to increment tallies with a delayed group filter. void -score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) +score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index, + std::vector& filter_matches) { // Save the original delayed group bin auto& tally {*model::tallies[i_tally]}; auto i_filt = tally.filters(tally.delayedgroup_filter_); - auto& dg_match {simulation::filter_matches[i_filt]}; + auto& dg_match {filter_matches[i_filt]}; auto i_bin = dg_match.i_bin_; auto original_bin = dg_match.bins_[i_bin]; dg_match.bins_[i_bin] = d_bin; @@ -158,7 +160,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) auto filter_index = 0; for (auto i = 0; i < tally.filters().size(); ++i) { auto i_filt = tally.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(i); } @@ -314,12 +316,12 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, //! neutrons produced with different energies. void -score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) +score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) { auto& tally {*model::tallies[i_tally]}; auto i_eout_filt = tally.filters()[tally.energyout_filter_]; - auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_; - auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin]; + auto i_bin = p->filter_matches_[i_eout_filt].i_bin_; + auto bin_energyout = p->filter_matches_[i_eout_filt].bins_[i_bin]; const EnergyoutFilter& eo_filt {*dynamic_cast(model::tally_filters[i_eout_filt].get())}; @@ -331,8 +333,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) // loop over number of particles banked for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; - const auto& bank = simulation::fission_bank[i_bank]; + const auto& bank = p->nu_bank_[i]; // get the delayed group auto g = bank.delayed_group; @@ -356,7 +357,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) g_out = eo_filt.n_bins() - g_out - 1; // change outgoing energy bin - simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out; + p->filter_matches_[i_eout_filt].bins_[i_bin] = g_out; } else { @@ -373,7 +374,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } else { auto i_match = lower_bound_index(eo_filt.bins().begin(), eo_filt.bins().end(), E_out); - simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match; + p->filter_matches_[i_eout_filt].bins_[i_bin] = i_match; } } @@ -387,7 +388,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) int filter_index = 0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); } @@ -415,13 +416,13 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_weight *= match.weights_[i_bin]; } score_fission_delayed_dg(i_tally, d_bin, score*filter_weight, - i_score); + i_score, p->filter_matches_); } } @@ -433,7 +434,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -447,7 +448,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } // Reset outgoing energy bin and score index - simulation::filter_matches[i_eout_filt].bins_[i_bin] = bin_energyout; + p->filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; } //! Update tally results for continuous-energy tallies with any estimator. @@ -777,7 +778,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, * p->neutron_xs_[p->event_nuclide_].fission / p->neutron_xs_[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } continue; } else { @@ -811,7 +812,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups()[d_bin]; score = simulation::keff * p->wgt_bank_ / p->n_bank_ * p->n_delayed_bank_[d-1] * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -836,7 +837,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->neutron_xs_[i_nuclide].fission * yield * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -866,7 +867,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->neutron_xs_[j_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } } } @@ -917,7 +918,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, / p->neutron_xs_[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } continue; } else { @@ -953,8 +954,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, // contribution to the fission bank to the score. score = 0.; for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; - const auto& bank = simulation::fission_bank[i_bank]; + const auto& bank = p->nu_bank_[i]; auto g = bank.delayed_group; if (g != 0) { const auto& nuc {*data::nuclides[p->event_nuclide_]}; @@ -971,7 +971,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups()[d_bin]; if (d == g) score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } score = 0.; } @@ -996,7 +996,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto rate = rxn.products_[d].decay_rate_; score = p->neutron_xs_[i_nuclide].fission * yield * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1037,7 +1037,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->neutron_xs_[j_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } } } @@ -1225,8 +1225,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, // We need to substract the energy of the secondary particles since // they will be transported individually later for (auto i = 0; i < p->n_bank_second_; ++i) { - auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i; - const auto& bank = simulation::secondary_bank[i_bank]; + auto i_bank = p->secondary_bank_.size() - p->n_bank_second_ + i; + const auto& bank = p->secondary_bank_[i_bank]; if (bank.particle == Particle::Type::photon || bank.particle == Particle::Type::neutron) { score -= bank.E; @@ -1342,7 +1342,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, //! argument is really just used for filter weights. void -score_general_mg(const Particle* p, int i_tally, int start_index, +score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { auto& tally {*model::tallies[i_tally]}; @@ -1723,7 +1723,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1764,7 +1764,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1795,7 +1795,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = flux * macro_xs.get_xs( MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1838,7 +1838,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1876,8 +1876,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // contribution to the fission bank to the score. score = 0.; for (auto i = 0; i < p->n_bank_; ++i) { - auto i_bank = simulation::fission_bank.size() - p->n_bank_ + i; - const auto& bank = simulation::fission_bank[i_bank]; + const auto& bank = p->nu_bank_[i]; auto d = bank.delayed_group - 1; if (d != -1) { if (i_nuclide >= 0) { @@ -1899,7 +1898,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, auto dg = filt.groups()[d_bin]; if (dg == d + 1) score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } score = 0.; } @@ -1929,7 +1928,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -2046,7 +2045,7 @@ void score_analog_tally_ce(Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2090,11 +2089,11 @@ void score_analog_tally_ce(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } -void score_analog_tally_mg(const Particle* p) +void score_analog_tally_mg(Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2103,7 +2102,7 @@ void score_analog_tally_mg(const Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2135,7 +2134,7 @@ void score_analog_tally_mg(const Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } @@ -2152,7 +2151,7 @@ score_tracklength_tally(Particle* p, double distance) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2199,7 +2198,7 @@ score_tracklength_tally(Particle* p, double distance) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } @@ -2220,7 +2219,7 @@ void score_collision_tally(Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2263,12 +2262,12 @@ void score_collision_tally(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } void -score_surface_tally(const Particle* p, const std::vector& tallies) +score_surface_tally(Particle* p, const std::vector& tallies) { // No collision, so no weight change when survival biasing double flux = p->wgt_; @@ -2280,7 +2279,7 @@ score_surface_tally(const Particle* p, const std::vector& tallies) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2307,7 +2306,7 @@ score_surface_tally(const Particle* p, const std::vector& tallies) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } diff --git a/src/timer.cpp b/src/timer.cpp index 2303c6090..616140aa9 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -16,6 +16,7 @@ Timer time_finalize; Timer time_inactive; Timer time_initialize; Timer time_read_xs; +Timer time_sample_source; Timer time_tallies; Timer time_total; Timer time_transport; @@ -68,6 +69,7 @@ void reset_timers() simulation::time_inactive.reset(); simulation::time_initialize.reset(); simulation::time_read_xs.reset(); + simulation::time_sample_source.reset(); simulation::time_tallies.reset(); simulation::time_total.reset(); simulation::time_transport.reset(); diff --git a/src/track_output.cpp b/src/track_output.cpp index 50d6d1ef6..5e071be01 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -13,37 +13,27 @@ #include #include -// Explicit vector template specialization of threadprivate variable outside of -// the openmc namespace for the picky Intel compiler. -template class std::vector>; - namespace openmc { //============================================================================== // Global variables //============================================================================== -// Forward declaration needed in order to declare tracks as threadprivate -extern std::vector> tracks; -#pragma omp threadprivate(tracks) - -std::vector> tracks; - //============================================================================== // Non-member functions //============================================================================== -void add_particle_track() +void add_particle_track(Particle& p) { - tracks.emplace_back(); + p.tracks_.emplace_back(); } -void write_particle_track(const Particle& p) +void write_particle_track(Particle& p) { - tracks.back().push_back(p.r()); + p.tracks_.back().push_back(p.r()); } -void finalize_particle_track(const Particle& p) +void finalize_particle_track(Particle& p) { std::stringstream filename; filename << settings::path_output << "track_" << simulation::current_batch @@ -51,7 +41,7 @@ void finalize_particle_track(const Particle& p) // Determine number of coordinates for each particle std::vector n_coords; - for (auto& coords : tracks) { + for (auto& coords : p.tracks_) { n_coords.push_back(coords.size()); } @@ -60,10 +50,10 @@ void finalize_particle_track(const Particle& p) hid_t file_id = file_open(filename.str().c_str(), 'w'); write_attribute(file_id, "filetype", "track"); write_attribute(file_id, "version", VERSION_TRACK); - write_attribute(file_id, "n_particles", tracks.size()); + write_attribute(file_id, "n_particles", p.tracks_.size()); write_attribute(file_id, "n_coords", n_coords); - for (int i = 1; i <= tracks.size(); ++i) { - const auto& t {tracks[i-1]}; + for (auto i = 1; i <= p.tracks_.size(); ++i) { + const auto& t {p.tracks_[i-1]}; size_t n = t.size(); xt::xtensor data({n,3}); for (int j = 0; j < n; ++j) { @@ -78,7 +68,7 @@ void finalize_particle_track(const Particle& p) } // Clear particle tracks - tracks.clear(); + p.tracks_.clear(); } } // namespace openmc diff --git a/tests/unit_tests/test_geometry.py b/tests/unit_tests/test_geometry.py index 70d90b615..64a97e0a1 100644 --- a/tests/unit_tests/test_geometry.py +++ b/tests/unit_tests/test_geometry.py @@ -301,3 +301,50 @@ def test_rotation_matrix(): assert geom.find((0.0, 0.5, 0.0))[-1] == c3 assert geom.find((0.0, -0.5, 0.0))[-1] == c1 assert geom.find((0.0, -1.5, 0.0))[-1] == c2 + +def test_remove_redundant_surfaces(): + """Test ability to remove redundant surfaces""" + + m1 = openmc.Material() + m1.add_nuclide('U235', 1.0, 'wo') + m1.add_nuclide('O16', 2.0, 'wo') + m1.set_density('g/cm3', 10.0) + + m2 = openmc.Material() + m2.add_element('Zr', 1.0) + m2.set_density('g/cm3', 2.0) + + m3 = openmc.Material() + m3.add_nuclide('H1', 2.0) + m3.add_nuclide('O16', 1.0) + m3.set_density('g/cm3', 1.0) + + def get_cyl_cell(r1, r2, z1, z2, fill): + """Create a finite height cylindrical cell with a fill""" + + cyl2 = openmc.ZCylinder(r=r2) + zplane1 = openmc.ZPlane(z1) + zplane2 = openmc.ZPlane(z2) + if np.isclose(r1, 0): + region = -cyl2 & +zplane1 & -zplane2 + else: + cyl1 = openmc.ZCylinder(r=r1) + region = +cyl1 & -cyl2 & +zplane1 & -zplane2 + return openmc.Cell(region=region, fill=fill) + + r1, r2, r3 = 1., 2., 3. + z1, z2 = -2., 2. + fuel = get_cyl_cell(0, r1, z1, z2, m1) + clad = get_cyl_cell(r1, r2, z1, z2, m2) + water = get_cyl_cell(r2, r3, z1, z2, m3) + root = openmc.Universe(cells=[fuel, clad, water]) + geom = openmc.Geometry(root) + + # There should be 6 redundant surfaces in this geometry + n_redundant_surfs = len(geom.get_redundant_surfaces().keys()) + assert n_redundant_surfs == 6 + # Remove redundant surfaces + geom.remove_redundant_surfaces() + # There should be 0 remaining redundant surfaces + n_redundant_surfs = len(geom.get_redundant_surfaces().keys()) + assert n_redundant_surfs == 0