diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 17f230866..b7fab9ac6 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -4,23 +4,13 @@ #include #include -namespace openmc { - -struct Bank { - double wgt; - double xyz[3]; - double uvw[3]; - double E; - int delayed_group; - int particle; -}; - -} // namespace openmc +#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; +extern template class std::vector; namespace openmc { @@ -32,10 +22,10 @@ namespace simulation { extern "C" int64_t n_bank; -extern std::vector source_bank; -extern std::vector fission_bank; +extern std::vector source_bank; +extern std::vector fission_bank; #ifdef _OPENMP -extern std::vector master_fission_bank; +extern std::vector master_fission_bank; #endif #pragma omp threadprivate(fission_bank, n_bank) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 599daf81c..9bfe4aa9e 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -6,23 +6,7 @@ #include #ifdef __cplusplus -#include "openmc/bank.h" - extern "C" { - int openmc_fission_bank(openmc::Bank** ptr, int64_t* n); - int openmc_source_bank(openmc::Bank** ptr, int64_t* n); -#else - struct Bank { - double wgt; - double xyz[3]; - double uvw[3]; - double E; - int delayed_group; - int particle; - }; - - int openmc_fission_bank(struct Bank** ptr, int64_t* n); - int openmc_source_bank(struct Bank** ptr, int64_t* n); #endif int openmc_calculate_volumes(); @@ -44,6 +28,7 @@ extern "C" { int openmc_filter_set_id(int32_t index, int32_t id); int openmc_finalize(); int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance); + int openmc_fission_bank(void** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); void openmc_get_filter_next_id(int32_t* id); @@ -90,6 +75,7 @@ extern "C" { void openmc_set_seed(int64_t new_seed); int openmc_simulation_finalize(); int openmc_simulation_init(); + int openmc_source_bank(void** ptr, int64_t* n); int openmc_spatial_legendre_filter_get_order(int32_t index, int* order); int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max); int openmc_spatial_legendre_filter_set_order(int32_t index, int order); diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 6343af605..f436592cd 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -3,6 +3,7 @@ #include #include +#include // for unique_ptr #include #include #include @@ -44,15 +45,11 @@ class Cell; class Universe; namespace model { + extern std::vector> cells; + extern std::unordered_map cell_map; -extern "C" int32_t n_cells; - -extern std::vector cells; -extern std::unordered_map cell_map; - -extern std::vector universes; -extern std::unordered_map universe_map; - + extern std::vector> universes; + extern std::unordered_map universe_map; } // namespace model //============================================================================== diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 077d789eb..e2662b71e 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -332,7 +332,17 @@ void read_dataset(hid_t obj_id, const char* name, xt::xtensor& arr, bool i // Copy into xtensor arr = xarr; +} +// overload for Position +inline void +read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false) +{ + std::array x; + read_dataset(obj_id, name, x, indep); + r.x = x[0]; + r.y = x[1]; + r.z = x[2]; } template diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 68728de5f..9c8b8af68 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -3,6 +3,7 @@ #include #include +#include // for unique_ptr #include #include #include @@ -33,10 +34,8 @@ enum class LatticeType { class Lattice; namespace model { - -extern std::vector lattices; -extern std::unordered_map lattice_map; - + extern std::vector> lattices; + extern std::unordered_map lattice_map; } // namespace model //============================================================================== @@ -105,7 +104,7 @@ public: //! \brief Find the lattice tile indices for a given point. //! \param r A 3D Cartesian coordinate. //! \return An array containing the indices of a lattice tile. - virtual std::array get_indices(Position r) const = 0; + virtual std::array get_indices(Position r, Direction u) const = 0; //! \brief Get coordinates local to a lattice tile. //! \param r A 3D Cartesian coordinate. @@ -213,7 +212,7 @@ public: std::pair> distance(Position r, Direction u, const std::array& i_xyz) const; - std::array get_indices(Position r) const; + std::array get_indices(Position r, Direction u) const; Position get_local_position(Position r, const std::array i_xyz) const; @@ -253,7 +252,7 @@ public: std::pair> distance(Position r, Direction u, const std::array& i_xyz) const; - std::array get_indices(Position r) const; + std::array get_indices(Position r, Direction u) const; Position get_local_position(Position r, const std::array i_xyz) const; diff --git a/include/openmc/material.h b/include/openmc/material.h index f43db5928..8a09dd9a2 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -23,7 +23,7 @@ class Material; namespace model { -extern std::vector materials; +extern std::vector> materials; extern std::unordered_map material_map; } // namespace model diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index e07374df3..39b703ab5 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -73,6 +73,8 @@ extern "C" double evaluate_legendre(int n, const double data[], double x); extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); +void calc_rn(int n, Direction u, double rn[]); + //============================================================================== //! Calculate the n-th order modified Zernike polynomial moment for a given //! angle (rho, theta) location on the unit disk. @@ -131,7 +133,7 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]); extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi); -Direction rotate_angle(Direction u, double mu, double* phi); +Direction rotate_angle(Direction u, double mu, const double* phi); //============================================================================== //! Samples an energy from the Maxwell fission distribution based on a direct diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index eed0d52f4..b13f9568a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -101,7 +101,7 @@ public: //! \param[in] energies Array of energies //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xarray count_sites(int64_t n, const Bank* bank, + xt::xarray count_sites(int64_t n, const Particle::Bank* bank, int n_energy, const double* energies, bool* outside) const; int id_ {-1}; //!< User-specified ID diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 145fbcf05..1b4062da9 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -158,12 +158,12 @@ class Mgxs { //! //! @param gin Incoming energy group. //! @param sqrtkT Temperature of the material. - //! @param uvw Incoming particle direction. + //! @param u Incoming particle direction. //! @param total_xs Resultant total cross section. //! @param abs_xs Resultant absorption cross section. //! @param nu_fiss_xs Resultant nu-fission cross section. void - calculate_xs(int gin, double sqrtkT, const double uvw[3], + calculate_xs(int gin, double sqrtkT, Direction u, double& total_xs, double& abs_xs, double& nu_fiss_xs); //! \brief Sets the temperature index in cache given a temperature @@ -174,9 +174,9 @@ class Mgxs { //! \brief Sets the angle index in cache given a direction //! - //! @param uvw Incoming particle direction. + //! @param u Incoming particle direction. void - set_angle_index(const double uvw[3]); + set_angle_index(Direction u); }; } // namespace openmc diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index fa65dd16f..bf72d675b 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -48,7 +48,7 @@ void read_mg_cross_sections_header(); //============================================================================== extern "C" void -calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], +calculate_xs_c(int i_mat, int gin, double sqrtkT, Direction u, double& total_xs, double& abs_xs, double& nu_fiss_xs); double diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 32d6c843b..2cb30a44f 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -9,7 +9,7 @@ #include #include -#include "openmc/capi.h" +#include "openmc/position.h" namespace openmc { @@ -33,20 +33,19 @@ constexpr int MAX_LOST_PARTICLES {10}; // Maximum number of lost particles, relative to the total number of particles constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; -//! Particle types -enum class ParticleType { - neutron, photon, electron, positron -}; +//============================================================================== +// Class declarations +//============================================================================== struct LocalCoord { + Position r; //!< particle position + Direction u; //!< particle direction int cell {-1}; int universe {-1}; int lattice {-1}; int lattice_x {-1}; int lattice_y {-1}; int lattice_z {-1}; - double xyz[3]; //!< particle position - double uvw[3]; //!< particle direction bool rotated {false}; //!< Is the level rotated? //! clear data from a single coordinate level @@ -57,70 +56,105 @@ struct LocalCoord { //! State of a particle being transported through geometry //============================================================================ -struct Particle { - int64_t id; //!< Unique ID - int type; //!< Particle type (n, p, e, etc.) +class Particle { +public: + //! Particle types + enum class Type { + neutron, photon, electron, positron + }; - int n_coord; //!< number of current coordinate levels - int cell_instance; //!< offset for distributed properties - LocalCoord coord[MAX_COORD]; //!< coordinates for all levels + //! Saved ("banked") state of a particle + struct Bank { + Position r; + Direction u; + double E; + double wgt; + int delayed_group; + Type particle; + }; + + // Constructors + Particle(); + + int64_t id_; //!< Unique ID + Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.) + + int n_coord_ {1}; //!< number of current coordinate levels + int cell_instance_; //!< offset for distributed properties + LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels // Particle coordinates before crossing a surface - int last_n_coord; //!< number of current coordinates - int last_cell[MAX_COORD]; //!< coordinates for all levels + int n_coord_last_ {1}; //!< number of current coordinates + int cell_last_[MAX_COORD]; //!< coordinates for all levels // Energy data - double E; //!< post-collision energy in eV - double last_E; //!< pre-collision energy in eV - int g; //!< post-collision energy group (MG only) - int last_g; //!< pre-collision energy group (MG only) + double E_; //!< post-collision energy in eV + double E_last_; //!< pre-collision energy in eV + int g_ {0}; //!< post-collision energy group (MG only) + int g_last_; //!< pre-collision energy group (MG only) // Other physical data - double wgt; //!< particle weight - double mu; //!< angle of scatter - bool alive; //!< is particle alive? + double wgt_ {1.0}; //!< particle weight + double mu_; //!< angle of scatter + bool alive_ {true}; //!< is particle alive? // Other physical data - double last_xyz_current[3]; //!< coordinates of the last collision or - //!< reflective/periodic surface crossing for - //!< current tallies - double last_xyz[3]; //!< previous coordinates - double last_uvw[3]; //!< previous direction coordinates - double last_wgt; //!< pre-collision particle weight - double absorb_wgt; //!< weight absorbed for survival biasing + Position r_last_current_; //!< coordinates of the last collision or + //!< reflective/periodic surface crossing for + //!< current tallies + Position r_last_; //!< previous coordinates + Direction u_last_; //!< previous direction coordinates + double wgt_last_ {1.0}; //!< pre-collision particle weight + double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing // What event took place - bool fission; //!< did particle cause implicit fission - int event; //!< scatter, absorption - int event_nuclide; //!< index in nuclides array - int event_MT; //!< reaction MT - int delayed_group; //!< delayed group + bool fission_ {false}; //!< did particle cause implicit fission + int event_; //!< scatter, absorption + int event_nuclide_; //!< index in nuclides array + int event_mt_; //!< reaction MT + int delayed_group_ {0}; //!< delayed group // Post-collision physical data - int n_bank; //!< number of fission sites banked - double wgt_bank; //!< weight of fission sites banked - int n_delayed_bank[MAX_DELAYED_GROUPS]; //!< number of delayed fission + int n_bank_ {0}; //!< number of fission sites banked + double wgt_bank_ {0.0}; //!< weight of fission sites banked + int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission //!< sites banked // Indices for various arrays - int surface; //!< index for surface particle is on - int cell_born; //!< index for cell particle was born in - int material; //!< index for current material - int last_material; //!< index for last material + 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 // Temperature of current cell - double sqrtkT; //!< sqrt(k_Boltzmann * temperature) in eV - double last_sqrtkT; //!< last temperature + double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV + double sqrtkT_last_ {0.0}; //!< last temperature // Statistical data - int n_collision; //!< number of collisions + int n_collision_ {0}; //!< number of collisions // Track output - bool write_track {false}; + bool write_track_ {false}; // Secondary particles created - int64_t n_secondary {}; - Bank secondary_bank[MAX_SECONDARY]; + int64_t n_secondary_ {}; + Bank secondary_bank_[MAX_SECONDARY]; + + // Accessors for position in global coordinates + Position& r() { return coord_[0].r; } + const Position& r() const { return coord_[0].r; } + + // Accessors for position in local coordinates + Position& r_local() { return coord_[n_coord_ - 1].r; } + const Position& r_local() const { return coord_[n_coord_ - 1].r; } + + // Accessors for direction in global coordinates + Direction& u() { return coord_[0].u; } + const Direction& u() const { return coord_[0].u; } + + // Accessors for direction in local coordinates + Direction& u_local() { return coord_[n_coord_ - 1].u; } + const Direction& u_local() const { return coord_[n_coord_ - 1].u; } //! resets all coordinate levels for the particle void clear(); @@ -129,14 +163,10 @@ struct Particle { // //! stores the current phase space attributes of the particle in the //! secondary bank and increments the number of sites in the secondary bank. - //! \param uvw Direction of the secondary particle + //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] //! \param type Particle type - //! \param run_CE Whether continuous-energy data is being used - void create_secondary(const double* uvw, double E, int type, bool run_CE); - - //! sets default attributes for a particle - void initialize(); + void create_secondary(Direction u, double E, Type type); //! initialize from a source site // diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 0c2b399d1..a7dd20aad 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -47,7 +47,7 @@ int sample_nuclide(const 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, - Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); + Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); int sample_element(Particle* p); @@ -60,18 +60,18 @@ void absorption(Particle* p, int i_nuclide); void scatter(Particle*, int i_nuclide); //! Treats the elastic scattering of a neutron with a target. -void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, - double* uvw, double* mu_lab, double* wgt); +void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, + Direction& u, double& mu_lab); -void sab_scatter(int i_nuclide, int i_sab, double* E, - double* uvw, double* mu); +void sab_scatter(int i_nuclide, int i_sab, double& E, + Direction& u, double& mu); //! samples the target velocity. The constant cross section free gas model is //! the default method. Methods for correctly accounting for the energy //! dependence of cross sections in treating resonance elastic scattering such //! as the DBRC and a new, accelerated scheme are also implemented here. Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, double* wgt); + Direction v_neut, double xs_eff, double kT); //! samples a target velocity based on the free gas scattering formulation, used //! by most Monte Carlo codes, in which cross section is assumed to be constant @@ -79,7 +79,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, //! FRA-TM-123. Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT); -void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site); +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site); //! handles all reactions with a single secondary neutron (other than fission), //! i.e. level scattering, (n,np), (n,na), etc. diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 511a0658b..cdf1f7089 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -35,7 +35,7 @@ scatter(Particle* p); //! \param size_bank Number of particles currently in the bank //! \param bank_array_size Allocated size of the bank void -create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, +create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_array_size); //! \brief Handles an absorption event diff --git a/include/openmc/position.h b/include/openmc/position.h index 0deb298dd..728763f25 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -1,6 +1,7 @@ #ifndef OPENMC_POSITION_H #define OPENMC_POSITION_H +#include #include // for sqrt #include // for out_of_range #include @@ -16,7 +17,8 @@ struct Position { Position() = default; Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; - Position(const std::vector xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + Position(const std::vector& xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + Position(const std::array& xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; // Unary operators Position& operator+=(Position); @@ -27,6 +29,7 @@ struct Position { Position& operator*=(double); Position& operator/=(Position); Position& operator/=(double); + Position operator-() const; const double& operator[](int i) const { switch (i) { diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index cd9f4636e..9377b5b6a 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -44,7 +44,7 @@ public: //! \param[out] mu Outgoing cosine with respect to current direction void sample(double E_in, double& E_out, double& mu) const; - ParticleType particle_; //!< Particle type + Particle::Type particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s] std::unique_ptr yield_; //!< Yield as a function of energy diff --git a/include/openmc/source.h b/include/openmc/source.h index 1fef129ff..bdd72e5bf 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -9,7 +9,6 @@ #include "pugixml.hpp" -#include "openmc/bank.h" #include "openmc/distribution_multi.h" #include "openmc/distribution_spatial.h" #include "openmc/particle.h" @@ -40,12 +39,12 @@ public: //! Sample from the external source distribution //! \return Sampled site - Bank sample() const; + Particle::Bank sample() const; // Properties double strength() const { return strength_; } private: - ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted + Particle::Type particle_ {Particle::Type::neutron}; //!< Type of particle emitted double strength_ {1.0}; //!< Source strength UPtrSpace space_; //!< Spatial distribution UPtrAngle angle_; //!< Angular distribution @@ -62,7 +61,7 @@ extern "C" void initialize_source(); //! Sample a site from all external source distributions in proportion to their //! source strength //! \return Sampled source site -Bank sample_external_source(); +Particle::Bank sample_external_source(); //! Fill source bank at end of generation for fixed source simulations void fill_source_bank_fixedsource(); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 69944ef92..e1acefa23 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -1,9 +1,10 @@ #ifndef OPENMC_SURFACE_H #define OPENMC_SURFACE_H -#include +#include // for unique_ptr #include // For numeric_limits #include +#include #include #include "hdf5.h" @@ -35,10 +36,8 @@ extern "C" const int BC_PERIODIC; class Surface; namespace model { - -extern std::vector surfaces; -extern std::map surface_map; - + extern std::vector> surfaces; + extern std::unordered_map surface_map; } // namespace model //============================================================================== diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h index 9f9328ace..268c19c6d 100644 --- a/include/openmc/tallies/filter_particle.h +++ b/include/openmc/tallies/filter_particle.h @@ -3,6 +3,7 @@ #include +#include "openmc/particle.h" #include "openmc/tallies/filter.h" namespace openmc { @@ -27,7 +28,7 @@ public: std::string text_label(int bin) const override; - std::vector particles_; + std::vector particles_; }; } // namespace openmc diff --git a/openmc/capi/core.py b/openmc/capi/core.py index dcca16454..922877c8e 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -13,10 +13,10 @@ import openmc.capi class _Bank(Structure): - _fields_ = [('wgt', c_double), - ('xyz', c_double*3), - ('uvw', c_double*3), + _fields_ = [('r', c_double*3), + ('u', c_double*3), ('E', c_double), + ('wgt', c_double), ('delayed_group', c_int)] diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 03bc22eae..d1ea77fa0 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1164,8 +1164,8 @@ class CMFDRun(object): energy = self._egrid ng = self._indices[3] - # Get xyz locations and energies of all particles in source bank - source_xyz = openmc.capi.source_bank()['xyz'] + # Get locations and energies of all particles in source bank + source_xyz = openmc.capi.source_bank()['r'] source_energies = openmc.capi.source_bank()['E'] # Convert xyz location to the CMFD mesh index @@ -1213,8 +1213,8 @@ class CMFDRun(object): outside = np.zeros(1, dtype=bool) count = np.zeros(self._sourcecounts.shape) - # Get xyz locations and energies of each particle in source bank - source_xyz = openmc.capi.source_bank()['xyz'] + # Get location and energy of each particle in source bank + source_xyz = openmc.capi.source_bank()['r'] source_energies = openmc.capi.source_bank()['E'] # Convert xyz location to mesh index and ravel index to scalar diff --git a/src/bank.cpp b/src/bank.cpp index 5f2f9cec5..8b7b43688 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -6,7 +6,7 @@ #include // Explicit template instantiation definition -template class std::vector; +template class std::vector; namespace openmc { @@ -18,10 +18,10 @@ namespace simulation { int64_t n_bank; -std::vector source_bank; -std::vector fission_bank; +std::vector source_bank; +std::vector fission_bank; #ifdef _OPENMP -std::vector master_fission_bank; +std::vector master_fission_bank; #endif } // namespace simulation @@ -46,7 +46,7 @@ void free_memory_bank() // C API //============================================================================== -extern "C" int openmc_source_bank(Bank** ptr, int64_t* n) +extern "C" int openmc_source_bank(void** ptr, int64_t* n) { if (simulation::source_bank.size() == 0) { set_errmsg("Source bank has not been allocated."); @@ -58,7 +58,7 @@ extern "C" int openmc_source_bank(Bank** ptr, int64_t* n) } } -extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n) +extern "C" int openmc_fission_bank(void** ptr, int64_t* n) { if (simulation::fission_bank.size() == 0) { set_errmsg("Fission bank has not been allocated."); diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index 2b44154da..3b9d8b591 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -28,20 +28,20 @@ std::vector ttb; void thick_target_bremsstrahlung(Particle& p, double* E_lost) { - if (p.material == MATERIAL_VOID) return; + if (p.material_ == MATERIAL_VOID) return; - int photon = static_cast(ParticleType::photon); - if (p.E < settings::energy_cutoff[photon]) return; + int photon = static_cast(Particle::Type::photon); + if (p.E_ < settings::energy_cutoff[photon]) return; // Get bremsstrahlung data for this material and particle type BremsstrahlungData* mat; - if (p.type == static_cast(ParticleType::positron)) { - mat = &model::materials[p.material]->ttb_->positron; + if (p.type_ == Particle::Type::positron) { + mat = &model::materials[p.material_]->ttb_->positron; } else { - mat = &model::materials[p.material]->ttb_->electron; + mat = &model::materials[p.material_]->ttb_->electron; } - double e = std::log(p.E); + double e = std::log(p.E_); auto n_e = data::ttb_e_grid.size(); // Find the lower bounding index of the incident electron energy @@ -108,8 +108,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon - int photon_ = static_cast(ParticleType::photon); - p.create_secondary(p.coord[0].uvw, w, photon_, true); + p.create_secondary(p.u(), w, Particle::Type::photon); *E_lost += w; } } diff --git a/src/cell.cpp b/src/cell.cpp index 2608d8f6f..d4712b9d1 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -22,13 +22,11 @@ namespace openmc { //============================================================================== namespace model { + std::vector> cells; + std::unordered_map cell_map; -std::vector cells; -std::unordered_map cell_map; - -std::vector universes; -std::unordered_map universe_map; - + std::vector> universes; + std::unordered_map universe_map; } // namespace model //============================================================================== @@ -653,7 +651,7 @@ void read_cells(pugi::xml_node node) // Loop over XML cell elements and populate the array. model::cells.reserve(n_cells); for (pugi::xml_node cell_node : node.children("cell")) { - model::cells.push_back(new CSGCell(cell_node)); + model::cells.push_back(std::make_unique(cell_node)); } // Fill the cell map. @@ -674,7 +672,7 @@ void read_cells(pugi::xml_node node) int32_t uid = model::cells[i]->universe_; auto it = model::universe_map.find(uid); if (it == model::universe_map.end()) { - model::universes.push_back(new Universe()); + model::universes.push_back(std::make_unique()); model::universes.back()->id_ = uid; model::universes.back()->cells_.push_back(i); model::universe_map[uid] = model::universes.size() - 1; @@ -823,7 +821,7 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end) if (index_start) *index_start = model::cells.size(); if (index_end) *index_end = model::cells.size() + n - 1; for (int32_t i = 0; i < n; i++) { - model::cells.push_back(new CSGCell()); + model::cells.push_back(std::make_unique()); } return 0; } diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 7cbb0789c..e91cfc875 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -228,7 +228,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // Determine if minimum/maximum energy for this nuclide is greater/less // than the previous if (data::nuclides[i_nuclide]->grid_.size() >= 1) { - int neutron = static_cast(ParticleType::neutron); + int neutron = static_cast(Particle::Type::neutron); data::energy_min[neutron] = std::max(data::energy_min[neutron], data::nuclides[i_nuclide]->grid_[0].energy.front()); data::energy_max[neutron] = std::min(data::energy_max[neutron], @@ -262,7 +262,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // the previous const auto& elem {data::elements.back()}; if (elem.energy_.size() >= 1) { - int photon = static_cast(ParticleType::photon); + int photon = static_cast(Particle::Type::photon); int n = elem.energy_.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(elem.energy_(1))); @@ -321,7 +321,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, for (auto& nuc : data::nuclides) { nuc->init_grid(); } - int neutron = static_cast(ParticleType::neutron); + int neutron = static_cast(Particle::Type::neutron); simulation::log_spacing = std::log(data::energy_max[neutron] / data::energy_min[neutron]) / settings::n_log_bins; @@ -329,7 +329,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // Determine if minimum/maximum energy for bremsstrahlung is greater/less // than the current minimum/maximum if (data::ttb_e_grid.size() >= 1) { - int photon = static_cast(ParticleType::photon); + int photon = static_cast(Particle::Type::photon); int n_e = data::ttb_e_grid.size(); data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); @@ -345,7 +345,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // grid has not been allocated if (nuc->grid_.size() > 0) { double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(ParticleType::neutron); + int neutron = static_cast(Particle::Type::neutron); if (max_E == data::energy_max[neutron]) { write_message("Maximum neutron transport energy: " + std::to_string(data::energy_max[neutron]) + " eV for " + diff --git a/src/dagmc.cpp b/src/dagmc.cpp index f417c12f5..34027ce7c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -143,14 +143,14 @@ void load_dagmc_geometry() c->universe_ = dagmc_univ_id; // set to zero for now c->fill_ = C_NONE; // no fill, single universe - model::cells.push_back(c); + model::cells.emplace_back(c); model::cell_map[c->id_] = i; // Populate the Universe vector and dict auto it = model::universe_map.find(dagmc_univ_id); if (it == model::universe_map.end()) { - model::universes.push_back(new Universe()); - model::universes.back()-> id_ = dagmc_univ_id; + model::universes.push_back(std::make_unique()); + model::universes.back()->id_ = dagmc_univ_id; model::universes.back()->cells_.push_back(i); model::universe_map[dagmc_univ_id] = model::universes.size() - 1; } else { @@ -255,7 +255,6 @@ void load_dagmc_geometry() // initialize surface objects int n_surfaces = model::DAG->num_entities(2); - model::surfaces.resize(n_surfaces); for (int i = 0; i < n_surfaces; i++) { moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1); @@ -303,7 +302,7 @@ void load_dagmc_geometry() } // add to global array and map - model::surfaces[i] = s; + model::surfaces.emplace_back(s); model::surface_map[s->id_] = s->id_; } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ccf70ee2e..c995a02d7 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -132,7 +132,7 @@ void synchronize_bank() // Allocate temporary source bank int64_t index_temp = 0; - std::vector temp_sites(3*simulation::work); + std::vector temp_sites(3*simulation::work); for (int64_t i = 0; i < simulation::n_bank; ++i) { // If there are less than n_particles particles banked, automatically add @@ -608,7 +608,7 @@ double ufs_get_weight(const Particle* p) auto& m = model::meshes[settings::index_ufs_mesh]; // Determine indices on ufs mesh for current location - int mesh_bin = m->get_bin({p->coord[0].xyz}); + int mesh_bin = m->get_bin(p->r()); if (mesh_bin < 0) { p->write_restart(); fatal_error("Source site outside UFS mesh!"); diff --git a/src/geometry.cpp b/src/geometry.cpp index fb081c313..3c1c11b0e 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -34,21 +34,21 @@ std::vector overlap_check_count; extern "C" bool check_cell_overlap(Particle* p) { - int n_coord = p->n_coord; + int n_coord = p->n_coord_; // Loop through each coordinate level for (int j = 0; j < n_coord; j++) { - Universe& univ = *model::universes[p->coord[j].universe]; + Universe& univ = *model::universes[p->coord_[j].universe]; int n = univ.cells_.size(); // Loop through each cell on this level for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; - if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) { - if (index_cell != p->coord[j].cell) { + if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) { + if (index_cell != p->coord_[j].cell) { std::stringstream err_msg; err_msg << "Overlapping cells detected: " << c.id_ << ", " - << model::cells[p->coord[j].cell]->id_ << " on universe " + << model::cells[p->coord_[j].cell]->id_ << " on universe " << univ.id_; fatal_error(err_msg); } @@ -74,36 +74,36 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) i_cell = *it; // Make sure the search cell is in the same universe. - int i_universe = p->coord[p->n_coord-1].universe; + int i_universe = p->coord_[p->n_coord_-1].universe; if (model::cells[i_cell]->universe_ != i_universe) continue; // Check if this cell contains the particle. - Position r {p->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; - auto surf = p->surface; + Position r {p->r_local()}; + Direction u {p->u_local()}; + auto surf = p->surface_; if (model::cells[i_cell]->contains(r, u, surf)) { - p->coord[p->n_coord-1].cell = i_cell; + p->coord_[p->n_coord_-1].cell = i_cell; found = true; break; } } } else { - int i_universe = p->coord[p->n_coord-1].universe; + int i_universe = p->coord_[p->n_coord_-1].universe; const auto& cells {model::universes[i_universe]->cells_}; for (auto it = cells.cbegin(); it != cells.cend(); it++) { i_cell = *it; // Make sure the search cell is in the same universe. - int i_universe = p->coord[p->n_coord-1].universe; + int i_universe = p->coord_[p->n_coord_-1].universe; if (model::cells[i_cell]->universe_ != i_universe) continue; // Check if this cell contains the particle. - Position r {p->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; - auto surf = p->surface; + Position r {p->r_local()}; + Direction u {p->u_local()}; + auto surf = p->surface_; if (model::cells[i_cell]->contains(r, u, surf)) { - p->coord[p->n_coord-1].cell = i_cell; + p->coord_[p->n_coord_-1].cell = i_cell; found = true; break; } @@ -126,37 +126,37 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) // Find the distribcell instance number. if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) { int offset = 0; - for (int i = 0; i < p->n_coord; i++) { - Cell& c_i {*model::cells[p->coord[i].cell]}; + for (int i = 0; i < p->n_coord_; i++) { + const auto& c_i {*model::cells[p->coord_[i].cell]}; if (c_i.type_ == FILL_UNIVERSE) { offset += c_i.offset_[c.distribcell_index_]; } else if (c_i.type_ == FILL_LATTICE) { - Lattice& lat {*model::lattices[p->coord[i+1].lattice-1]}; - int i_xyz[3] {p->coord[i+1].lattice_x, - p->coord[i+1].lattice_y, - p->coord[i+1].lattice_z}; + auto& lat {*model::lattices[p->coord_[i+1].lattice]}; + int i_xyz[3] {p->coord_[i+1].lattice_x, + p->coord_[i+1].lattice_y, + p->coord_[i+1].lattice_z}; if (lat.are_valid_indices(i_xyz)) { offset += lat.offset(c.distribcell_index_, i_xyz); } } } - p->cell_instance = offset; + p->cell_instance_ = offset; } else { - p->cell_instance = 0; + p->cell_instance_ = 0; } // Set the material and temperature. - p->last_material = p->material; + p->material_last_ = p->material_; if (c.material_.size() > 1) { - p->material = c.material_[p->cell_instance]; + p->material_ = c.material_[p->cell_instance_]; } else { - p->material = c.material_[0]; + p->material_ = c.material_[0]; } - p->last_sqrtkT = p->sqrtkT; + p->sqrtkT_last_ = p->sqrtkT_; if (c.sqrtkT_.size() > 1) { - p->sqrtkT = c.sqrtkT_[p->cell_instance]; + p->sqrtkT_ = c.sqrtkT_[p->cell_instance_]; } else { - p->sqrtkT = c.sqrtkT_[0]; + p->sqrtkT_ = c.sqrtkT_[0]; } return true; @@ -166,44 +166,36 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) //! Found a lower universe, update this coord level then search the next. // Set the lower coordinate level universe. - p->coord[p->n_coord].universe = c.fill_; + p->coord_[p->n_coord_].universe = c.fill_; // Set the position and direction. - for (int i = 0; i < 3; i++) { - p->coord[p->n_coord].xyz[i] = p->coord[p->n_coord-1].xyz[i]; - p->coord[p->n_coord].uvw[i] = p->coord[p->n_coord-1].uvw[i]; - } + p->coord_[p->n_coord_].r = p->coord_[p->n_coord_-1].r; + p->coord_[p->n_coord_].u = p->coord_[p->n_coord_-1].u; // Apply translation. - p->coord[p->n_coord].xyz[0] -= c.translation_.x; - p->coord[p->n_coord].xyz[1] -= c.translation_.y; - p->coord[p->n_coord].xyz[2] -= c.translation_.z; + p->coord_[p->n_coord_].r -= c.translation_; // Apply rotation. if (!c.rotation_.empty()) { - auto x = p->coord[p->n_coord].xyz[0]; - auto y = p->coord[p->n_coord].xyz[1]; - auto z = p->coord[p->n_coord].xyz[2]; - p->coord[p->n_coord].xyz[0] = x*c.rotation_[3] + y*c.rotation_[4] - + z*c.rotation_[5]; - p->coord[p->n_coord].xyz[1] = x*c.rotation_[6] + y*c.rotation_[7] - + z*c.rotation_[8]; - p->coord[p->n_coord].xyz[2] = x*c.rotation_[9] + y*c.rotation_[10] - + z*c.rotation_[11]; - auto u = p->coord[p->n_coord].uvw[0]; - auto v = p->coord[p->n_coord].uvw[1]; - auto w = p->coord[p->n_coord].uvw[2]; - p->coord[p->n_coord].uvw[0] = u*c.rotation_[3] + v*c.rotation_[4] - + w*c.rotation_[5]; - p->coord[p->n_coord].uvw[1] = u*c.rotation_[6] + v*c.rotation_[7] - + w*c.rotation_[8]; - p->coord[p->n_coord].uvw[2] = u*c.rotation_[9] + v*c.rotation_[10] - + w*c.rotation_[11]; - p->coord[p->n_coord].rotated = true; + Position r = p->coord_[p->n_coord_].r; + p->coord_[p->n_coord_].r.x = r.x*c.rotation_[3] + r.y*c.rotation_[4] + + r.z*c.rotation_[5]; + p->coord_[p->n_coord_].r.y = r.x*c.rotation_[6] + r.y*c.rotation_[7] + + r.z*c.rotation_[8]; + p->coord_[p->n_coord_].r.z = r.x*c.rotation_[9] + r.y*c.rotation_[10] + + r.z*c.rotation_[11]; + Direction u = p->coord_[p->n_coord_].u; + p->coord_[p->n_coord_].u.x = u.x*c.rotation_[3] + u.y*c.rotation_[4] + + u.z*c.rotation_[5]; + p->coord_[p->n_coord_].u.y = u.x*c.rotation_[6] + u.y*c.rotation_[7] + + u.z*c.rotation_[8]; + p->coord_[p->n_coord_].u.z = u.x*c.rotation_[9] + u.y*c.rotation_[10] + + u.z*c.rotation_[11]; + p->coord_[p->n_coord_].rotated = true; } // Update the coordinate level and recurse. - ++p->n_coord; + ++p->n_coord_; return find_cell_inner(p, nullptr); } else if (c.type_ == FILL_LATTICE) { @@ -213,35 +205,28 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) Lattice& lat {*model::lattices[c.fill_]}; // Determine lattice indices. - Position r {p->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; - r += TINY_BIT * u; - auto i_xyz = lat.get_indices(r); + auto i_xyz = lat.get_indices(p->r_local(), p->u_local()); // Store lower level coordinates. - r = lat.get_local_position(p->coord[p->n_coord-1].xyz, i_xyz); - p->coord[p->n_coord].xyz[0] = r.x; - p->coord[p->n_coord].xyz[1] = r.y; - p->coord[p->n_coord].xyz[2] = r.z; - p->coord[p->n_coord].uvw[0] = u.x; - p->coord[p->n_coord].uvw[1] = u.y; - p->coord[p->n_coord].uvw[2] = u.z; + Position r = lat.get_local_position(p->r_local(), i_xyz); + p->coord_[p->n_coord_].r = r; + p->coord_[p->n_coord_].u = p->u_local(); // Set lattice indices. - p->coord[p->n_coord].lattice = c.fill_ + 1; - p->coord[p->n_coord].lattice_x = i_xyz[0]; - p->coord[p->n_coord].lattice_y = i_xyz[1]; - p->coord[p->n_coord].lattice_z = i_xyz[2]; + p->coord_[p->n_coord_].lattice = c.fill_; + p->coord_[p->n_coord_].lattice_x = i_xyz[0]; + p->coord_[p->n_coord_].lattice_y = i_xyz[1]; + p->coord_[p->n_coord_].lattice_z = i_xyz[2]; // Set the lower coordinate level universe. if (lat.are_valid_indices(i_xyz)) { - p->coord[p->n_coord].universe = lat[i_xyz]; + p->coord_[p->n_coord_].universe = lat[i_xyz]; } else { if (lat.outer_ != NO_OUTER_UNIVERSE) { - p->coord[p->n_coord].universe = lat.outer_; + p->coord_[p->n_coord_].universe = lat.outer_; } else { std::stringstream err_msg; - err_msg << "Particle " << p->id << " is outside lattice " + err_msg << "Particle " << p->id_ << " is outside lattice " << lat.id_ << " but the lattice has no defined outer " "universe."; warning(err_msg); @@ -250,7 +235,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } // Update the coordinate level and recurse. - ++p->n_coord; + ++p->n_coord_; return find_cell_inner(p, nullptr); } } @@ -264,22 +249,22 @@ extern "C" bool find_cell(Particle* p, bool use_neighbor_lists) { // Determine universe (if not yet set, use root universe). - int i_universe = p->coord[p->n_coord-1].universe; + int i_universe = p->coord_[p->n_coord_-1].universe; if (i_universe == C_NONE) { - p->coord[0].universe = model::root_universe; - p->n_coord = 1; + p->coord_[0].universe = model::root_universe; + p->n_coord_ = 1; i_universe = model::root_universe; } // Reset all the deeper coordinate levels. - for (int i = p->n_coord; i < MAX_COORD; i++) { - p->coord[i].reset(); + for (int i = p->n_coord_; i < MAX_COORD; i++) { + p->coord_[i].reset(); } if (use_neighbor_lists) { // Get the cell this particle was in previously. - auto coord_lvl = p->n_coord - 1; - auto i_cell = p->coord[coord_lvl].cell; + auto coord_lvl = p->n_coord_ - 1; + auto i_cell = p->coord_[coord_lvl].cell; Cell& c {*model::cells[i_cell]}; // Search for the particle in that cell's neighbor list. Return if we @@ -291,7 +276,7 @@ find_cell(Particle* p, bool use_neighbor_lists) // cells in this universe, and update the neighbor list if we find a new // neighboring cell. found = find_cell_inner(p, nullptr); - if (found) c.neighbors_.push_back(p->coord[coord_lvl].cell); + if (found) c.neighbors_.push_back(p->coord_[coord_lvl].cell); return found; } else { @@ -305,55 +290,52 @@ find_cell(Particle* p, bool use_neighbor_lists) extern "C" void cross_lattice(Particle* p, int lattice_translation[3]) { - Lattice& lat {*model::lattices[p->coord[p->n_coord-1].lattice-1]}; + auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]}; if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; msg << " Crossing lattice " << lat.id_ << ". Current position (" - << p->coord[p->n_coord-1].lattice_x << "," - << p->coord[p->n_coord-1].lattice_y << "," - << p->coord[p->n_coord-1].lattice_z << ")"; + << p->coord_[p->n_coord_-1].lattice_x << "," + << p->coord_[p->n_coord_-1].lattice_y << "," + << p->coord_[p->n_coord_-1].lattice_z << ")"; write_message(msg, 1); } // Set the lattice indices. - p->coord[p->n_coord-1].lattice_x += lattice_translation[0]; - p->coord[p->n_coord-1].lattice_y += lattice_translation[1]; - p->coord[p->n_coord-1].lattice_z += lattice_translation[2]; - std::array i_xyz {p->coord[p->n_coord-1].lattice_x, - p->coord[p->n_coord-1].lattice_y, - p->coord[p->n_coord-1].lattice_z}; + p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0]; + p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1]; + p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2]; + std::array i_xyz {p->coord_[p->n_coord_-1].lattice_x, + p->coord_[p->n_coord_-1].lattice_y, + p->coord_[p->n_coord_-1].lattice_z}; // Set the new coordinate position. - auto r = lat.get_local_position(p->coord[p->n_coord-2].xyz, i_xyz); - p->coord[p->n_coord-1].xyz[0] = r.x; - p->coord[p->n_coord-1].xyz[1] = r.y; - p->coord[p->n_coord-1].xyz[2] = r.z; + p->r_local() = lat.get_local_position(p->coord_[p->n_coord_-2].r, i_xyz); if (!lat.are_valid_indices(i_xyz)) { // The particle is outside the lattice. Search for it from the base coords. - p->n_coord = 1; + p->n_coord_ = 1; bool found = find_cell(p, 0); - if (!found && p->alive) { + if (!found && p->alive_) { std::stringstream err_msg; - err_msg << "Could not locate particle " << p->id + err_msg << "Could not locate particle " << p->id_ << " after crossing a lattice boundary"; p->mark_as_lost(err_msg); } } else { // Find cell in next lattice element. - p->coord[p->n_coord-1].universe = lat[i_xyz]; + p->coord_[p->n_coord_-1].universe = lat[i_xyz]; bool found = find_cell(p, 0); if (!found) { // A particle crossing the corner of a lattice tile may not be found. In // this case, search for it from the base coords. - p->n_coord = 1; + p->n_coord_ = 1; bool found = find_cell(p, 0); - if (!found && p->alive) { + if (!found && p->alive_) { std::stringstream err_msg; - err_msg << "Could not locate particle " << p->id + err_msg << "Could not locate particle " << p->id_ << " after crossing a lattice boundary"; p->mark_as_lost(err_msg); } @@ -377,21 +359,21 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, std::array level_lat_trans; // Loop over each coordinate level. - for (int i = 0; i < p->n_coord; i++) { - Position r {p->coord[i].xyz}; - Direction u {p->coord[i].uvw}; - Cell& c {*model::cells[p->coord[i].cell]}; + for (int i = 0; i < p->n_coord_; i++) { + Position r {p->coord_[i].r}; + Direction u {p->coord_[i].u}; + Cell& c {*model::cells[p->coord_[i].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_); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; // Find the distance to the next lattice tile crossing. - if (p->coord[i].lattice != F90_NONE) { - Lattice& lat {*model::lattices[p->coord[i].lattice-1]}; - std::array i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y, - p->coord[i].lattice_z}; + if (p->coord_[i].lattice != C_NONE) { + auto& lat {*model::lattices[p->coord_[i].lattice]}; + std::array i_xyz {p->coord_[i].lattice_x, p->coord_[i].lattice_y, + p->coord_[i].lattice_z}; //TODO: refactor so both lattice use the same position argument (which //also means the lat.type attribute can be removed) std::pair> lattice_distance; @@ -400,8 +382,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, lattice_distance = lat.distance(r, u, i_xyz); break; case LatticeType::hex: - Position r_hex {p->coord[i-1].xyz[0], p->coord[i-1].xyz[1], - p->coord[i].xyz[2]}; + Position r_hex {p->coord_[i-1].r.x, p->coord_[i-1].r.y, + p->coord_[i].r.z}; lattice_distance = lat.distance(r_hex, u, i_xyz); break; } @@ -410,7 +392,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, if (d_lat < 0) { std::stringstream err_msg; - err_msg << "Particle " << p->id + err_msg << "Particle " << p->id_ << " had a negative distance to a lattice boundary"; p->mark_as_lost(err_msg); } @@ -466,23 +448,20 @@ extern "C" int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) { Particle p; - p.initialize(); - std::copy(xyz, xyz + 3, p.coord[0].xyz); - p.coord[0].uvw[0] = 0.0; - p.coord[0].uvw[1] = 0.0; - p.coord[0].uvw[2] = 1.0; + p.r() = Position{xyz}; + p.u() = {0.0, 0.0, 1.0}; if (!find_cell(&p, false)) { std::stringstream msg; - msg << "Could not find cell at position (" << xyz[0] << ", " << xyz[1] - << ", " << xyz[2] << ")."; + msg << "Could not find cell at position (" << p.r().x << ", " << p.r().y + << ", " << p.r().z << ")."; set_errmsg(msg); return OPENMC_E_GEOMETRY; } - *index = p.coord[p.n_coord-1].cell; - *instance = p.cell_instance; + *index = p.coord_[p.n_coord_-1].cell; + *instance = p.cell_instance_; return 0; } diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 9cdc93752..c560743f8 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -66,7 +66,7 @@ void adjust_indices() { // Adjust material/fill idices. - for (Cell* c : model::cells) { + for (auto& c : model::cells) { if (c->fill_ != C_NONE) { int32_t id = c->fill_; auto search_univ = model::universe_map.find(id); @@ -102,7 +102,7 @@ adjust_indices() } // Change cell.universe values from IDs to indices. - for (Cell* c : model::cells) { + for (auto& c : model::cells) { auto search = model::universe_map.find(c->universe_); if (search != model::universe_map.end()) { c->universe_ = search->second; @@ -115,7 +115,7 @@ adjust_indices() } // Change all lattice universe values from IDs to indices. - for (Lattice* l : model::lattices) { + for (auto& l : model::lattices) { l->adjust_indices(); } } @@ -125,7 +125,7 @@ adjust_indices() void assign_temperatures() { - for (Cell* c : model::cells) { + for (auto& c : model::cells) { // Ignore non-material cells and cells with defined temperature. if (c->material_.size() == 0) continue; if (c->sqrtkT_.size() > 0) continue; @@ -234,12 +234,12 @@ find_root_universe() { // Find all the universes listed as a cell fill. std::unordered_set fill_univ_ids; - for (Cell* c : model::cells) { + for (const auto& c : model::cells) { fill_univ_ids.insert(c->fill_); } // Find all the universes contained in a lattice. - for (Lattice* lat : model::lattices) { + for (const auto& lat : model::lattices) { for (auto it = lat->begin(); it != lat->end(); ++it) { fill_univ_ids.insert(*it); } @@ -318,7 +318,7 @@ prepare_distribcell() // unique distribcell array index. int distribcell_index = 0; std::vector target_univ_ids; - for (Universe* u : model::universes) { + for (const auto& u : model::universes) { for (auto cell_indx : u->cells_) { if (distribcells.find(cell_indx) != distribcells.end()) { model::cells[cell_indx]->distribcell_index_ = distribcell_index; @@ -330,19 +330,19 @@ prepare_distribcell() // Allocate the cell and lattice offset tables. int n_maps = target_univ_ids.size(); - for (Cell* c : model::cells) { + for (auto& c : model::cells) { if (c->type_ != FILL_MATERIAL) { c->offset_.resize(n_maps, C_NONE); } } - for (Lattice* lat : model::lattices) { + for (auto& lat : model::lattices) { lat->allocate_offset_table(n_maps); } // Fill the cell and lattice offset tables. for (int map = 0; map < target_univ_ids.size(); map++) { auto target_univ_id = target_univ_ids[map]; - for (Universe* univ : model::universes) { + for (const auto& univ : model::universes) { int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. for (int32_t cell_indx : univ->cells_) { Cell& c = *model::cells[cell_indx]; @@ -525,15 +525,12 @@ maximum_levels(int32_t univ) void free_memory_geometry() { - for (Cell* c : model::cells) {delete c;} model::cells.clear(); model::cell_map.clear(); - for (Universe* u : model::universes) {delete u;} model::universes.clear(); model::universe_map.clear(); - for (Lattice* lat : model::lattices) {delete lat;} model::lattices.clear(); model::lattice_map.clear(); diff --git a/src/initialize.cpp b/src/initialize.cpp index 4bb4817ea..3aa31184f 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -100,17 +100,17 @@ void initialize_mpi(MPI_Comm intracomm) mpi::master = (mpi::rank == 0); // Create bank datatype - Bank b; + Particle::Bank b; MPI_Aint disp[6]; - MPI_Get_address(&b.wgt, &disp[0]); - MPI_Get_address(&b.xyz, &disp[1]); - MPI_Get_address(&b.uvw, &disp[2]); - MPI_Get_address(&b.E, &disp[3]); + 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]; - int blocks[] {1, 3, 3, 1, 1, 1}; + 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); MPI_Type_commit(&mpi::bank); diff --git a/src/lattice.cpp b/src/lattice.cpp index d8185de62..6e0f2331f 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -19,10 +19,8 @@ namespace openmc { //============================================================================== namespace model { - -std::vector lattices; -std::unordered_map lattice_map; - + std::vector> lattices; + std::unordered_map lattice_map; } //============================================================================== @@ -287,15 +285,38 @@ const //============================================================================== std::array -RectLattice::get_indices(Position r) const +RectLattice::get_indices(Position r, Direction u) const { - int ix {static_cast(std::ceil((r.x - lower_left_.x) / pitch_.x))-1}; - int iy {static_cast(std::ceil((r.y - lower_left_.y) / pitch_.y))-1}; - int iz; - if (is_3d_) { - iz = static_cast(std::ceil((r.z - lower_left_.z) / pitch_.z))-1; + // Determine x index, accounting for coincidence + double ix_ {(r.x - lower_left_.x) / pitch_.x}; + long ix_close {std::lround(ix_)}; + int ix; + if (std::abs(ix_ - ix_close) < FP_COINCIDENT) { + ix = (u.x > 0) ? ix_close : ix_close - 1; } else { - iz = 0; + ix = std::floor(ix_); + } + + // Determine y index, accounting for coincidence + double iy_ {(r.y - lower_left_.y) / pitch_.y}; + long iy_close {std::lround(iy_)}; + int iy; + if (std::abs(iy_ - iy_close) < FP_COINCIDENT) { + iy = (u.y > 0) ? iy_close : iy_close - 1; + } else { + iy = std::floor(iy_); + } + + // Determine z index, accounting for coincidence + int iz = 0; + if (is_3d_) { + double iz_ {(r.z - lower_left_.z) / pitch_.z}; + long iz_close {std::lround(iz_)}; + if (std::abs(iz_ - iz_close) < FP_COINCIDENT) { + iz = (u.z > 0) ? iz_close : iz_close - 1; + } else { + iz = std::floor(iz_); + } } return {ix, iy, iz}; } @@ -687,8 +708,13 @@ const //============================================================================== std::array -HexLattice::get_indices(Position r) const +HexLattice::get_indices(Position r, Direction u) const { + // The implementation for HexLattice currently doesn't use direction + // information. As a result, we move the position slightly forward to + // determine what lattice index the particle is most likely to be in. + r += TINY_BIT * u; + // Offset the xyz by the lattice center. Position r_o {r.x - center_.x, r.y - center_.y, r.z}; if (is_3d_) {r_o.z -= center_.z;} @@ -696,7 +722,7 @@ HexLattice::get_indices(Position r) const // Index the z direction. std::array out; if (is_3d_) { - out[2] = static_cast(std::ceil(r_o.z / pitch_[1] + 0.5 * n_axial_))-1; + out[2] = std::floor(r_o.z / pitch_[1] + 0.5 * n_axial_); } else { out[2] = 0; } @@ -704,9 +730,8 @@ HexLattice::get_indices(Position r) const // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); - out[0] = static_cast(std::floor(r_o.x - / (0.5*std::sqrt(3.0) * pitch_[0]))); - out[1] = static_cast(std::floor(alpha / pitch_[0])); + out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0])); + out[1] = std::floor(alpha / pitch_[0]); // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but // the array is offset so that the indices never go below 0). @@ -739,12 +764,12 @@ HexLattice::get_indices(Position r) const // Select the minimum squared distance which corresponds to the cell the // coordinates are in. if (k_min == 2) { - out[0] += 1; + ++out[0]; } else if (k_min == 3) { - out[1] += 1; + ++out[1]; } else if (k_min == 4) { - out[0] += 1; - out[1] += 1; + ++out[0]; + ++out[1]; } return out; @@ -867,10 +892,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const void read_lattices(pugi::xml_node node) { for (pugi::xml_node lat_node : node.children("lattice")) { - model::lattices.push_back(new RectLattice(lat_node)); + model::lattices.push_back(std::make_unique(lat_node)); } for (pugi::xml_node lat_node : node.children("hex_lattice")) { - model::lattices.push_back(new HexLattice(lat_node)); + model::lattices.push_back(std::make_unique(lat_node)); } // Fill the lattice map. diff --git a/src/material.cpp b/src/material.cpp index 8236a4ef0..e6b2ca458 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -37,7 +37,7 @@ namespace openmc { namespace model { -std::vector materials; +std::vector> materials; std::unordered_map material_map; } // namespace model @@ -652,19 +652,18 @@ void Material::calculate_xs(const Particle& p) const simulation::material_xs.fission = 0.0; simulation::material_xs.nu_fission = 0.0; - if (p.type == static_cast(ParticleType::neutron)) { + if (p.type_ == Particle::Type::neutron) { this->calculate_neutron_xs(p); - } else if (p.type == static_cast(ParticleType::photon)) { + } else if (p.type_ == Particle::Type::photon) { this->calculate_photon_xs(p); } } void Material::calculate_neutron_xs(const Particle& p) const { - int neutron = static_cast(ParticleType::neutron); - // Find energy index on energy grid - int i_grid = std::log(p.E/data::energy_min[neutron])/simulation::log_spacing; + int neutron = static_cast(Particle::Type::neutron); + int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing; // Determine if this material has S(a,b) tables bool check_sab = (thermal_tables_.size() > 0); @@ -691,7 +690,7 @@ void Material::calculate_neutron_xs(const Particle& p) const // If particle energy is greater than the highest energy for the // S(a,b) table, then don't use the S(a,b) table - if (p.E > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; + if (p.E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; // Increment position in thermal_tables_ ++j; @@ -709,12 +708,12 @@ void Material::calculate_neutron_xs(const Particle& p) const // Calculate microscopic cross section for this nuclide const auto& micro {simulation::micro_xs[i_nuclide]}; - if (p.E != micro.last_E - || p.sqrtkT != micro.last_sqrtkT + if (p.E_ != micro.last_E + || p.sqrtkT_ != micro.last_sqrtkT || i_sab != micro.index_sab || sab_frac != micro.sab_frac) { - data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E, i_grid, - p.sqrtkT, sab_frac); + data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E_, i_grid, + p.sqrtkT_, sab_frac); } // ====================================================================== @@ -748,8 +747,8 @@ void Material::calculate_photon_xs(const Particle& p) const // Calculate microscopic cross section for this nuclide const auto& micro {simulation::micro_photon_xs[i_element]}; - if (p.E != micro.last_E) { - data::elements[i_element].calculate_xs(p.E); + if (p.E_ != micro.last_E) { + data::elements[i_element].calculate_xs(p.E_); } // ======================================================================== @@ -895,7 +894,7 @@ void read_materials_xml() // Loop over XML material elements and populate the array. pugi::xml_node root = doc.document_element(); for (pugi::xml_node material_node : root.children("material")) { - model::materials.push_back(new Material(material_node)); + model::materials.push_back(std::make_unique(material_node)); } model::materials.shrink_to_fit(); @@ -915,7 +914,6 @@ void read_materials_xml() void free_memory_material() { - for (Material *mat : model::materials) {delete mat;} model::materials.clear(); model::material_map.clear(); } @@ -942,7 +940,7 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density) { int err = 0; if (index >= 0 && index < model::materials.size()) { - Material* m = model::materials[index]; + auto& m = model::materials[index]; // Check if nuclide is already in material for (int i = 0; i < m->nuclide_.size(); ++i) { @@ -1032,7 +1030,7 @@ extern "C" int openmc_material_get_volume(int32_t index, double* volume) { if (index >= 0 && index < model::materials.size()) { - Material* m = model::materials[index]; + auto& m = model::materials[index]; if (m->volume_ >= 0.0) { *volume = m->volume_; return 0; @@ -1112,7 +1110,7 @@ extern "C" int openmc_material_set_volume(int32_t index, double volume) { if (index >= 0 && index < model::materials.size()) { - Material* m = model::materials[index]; + auto& m {model::materials[index]}; if (volume >= 0.0) { m->volume_ = volume; return 0; @@ -1132,7 +1130,7 @@ openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end) if (index_start) *index_start = model::materials.size(); if (index_end) *index_end = model::materials.size() + n - 1; for (int32_t i = 0; i < n; i++) { - model::materials.push_back(new Material()); + model::materials.push_back(std::make_unique()); } return 0; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 919c305c4..54d0e9344 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -8,7 +8,8 @@ namespace openmc { // Mathematical methods //============================================================================== -double normal_percentile(double p) { +double normal_percentile(double p) +{ constexpr double p_low = 0.02425; constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2, 1.383577518672690e2, @@ -60,7 +61,8 @@ double normal_percentile(double p) { } -double t_percentile(double p, int df){ +double t_percentile(double p, int df) +{ double t; if (df == 1) { @@ -92,7 +94,8 @@ double t_percentile(double p, int df){ } -void calc_pn_c(int n, double x, double pnx[]) { +void calc_pn_c(int n, double x, double pnx[]) +{ pnx[0] = 1.; if (n >= 1) { pnx[1] = x; @@ -105,7 +108,8 @@ void calc_pn_c(int n, double x, double pnx[]) { } -double evaluate_legendre(int n, const double data[], double x) { +double evaluate_legendre(int n, const double data[], double x) +{ double pnx[n + 1]; double val = 0.0; calc_pn_c(n, x, pnx); @@ -116,16 +120,24 @@ double evaluate_legendre(int n, const double data[], double x) { } -void calc_rn_c(int n, const double uvw[3], double rn[]){ +void calc_rn_c(int n, const double uvw[3], double rn[]) +{ + Direction u {uvw}; + calc_rn(n, u, rn); +} + + +void calc_rn(int n, Direction u, double rn[]) +{ // rn[] is assumed to have already been allocated to the correct size // Store the cosine of the polar angle and the azimuthal angle - double w = uvw[2]; + double w = u.z; double phi; - if (uvw[0] == 0.) { + if (u.x == 0.) { phi = 0.; } else { - phi = std::atan2(uvw[1], uvw[0]); + phi = std::atan2(u.y, u.x); } // Store the shorthand of 1-w * w @@ -618,48 +630,44 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) { void rotate_angle_c(double uvw[3], double mu, const double* phi) { - // Copy original directional cosines - double u0 = uvw[0]; // original cosine in x direction - double v0 = uvw[1]; // original cosine in y direction - double w0 = uvw[2]; // original cosine in z direction + Direction u = rotate_angle({uvw}, mu, phi); + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; +} + +Direction rotate_angle(Direction u, double mu, const double* phi) +{ // Sample azimuthal angle in [0,2pi) if none provided double phi_; if (phi != nullptr) { phi_ = (*phi); } else { - phi_ = 2. * PI * prn(); + phi_ = 2.0*PI*prn(); } // Precompute factors to save flops double sinphi = std::sin(phi_); double cosphi = std::cos(phi_); - double a = std::sqrt(std::fmax(0., 1. - mu * mu)); - double b = std::sqrt(std::fmax(0., 1. - w0 * w0)); + double a = std::sqrt(std::fmax(0., 1. - mu*mu)); + double b = std::sqrt(std::fmax(0., 1. - u.z*u.z)); // Need to treat special case where sqrt(1 - w**2) is close to zero by // expanding about the v component rather than the w component if (b > 1e-10) { - uvw[0] = mu * u0 + a * (u0 * w0 * cosphi - v0 * sinphi) / b; - uvw[1] = mu * v0 + a * (v0 * w0 * cosphi + u0 * sinphi) / b; - uvw[2] = mu * w0 - a * b * cosphi; + return {mu*u.x + a*(u.x*u.z*cosphi - u.y*sinphi) / b, + mu*u.y + a*(u.y*u.z*cosphi + u.x*sinphi) / b, + mu*u.z - a*b*cosphi}; } else { - b = std::sqrt(1. - v0 * v0); - uvw[0] = mu * u0 + a * (u0 * v0 * cosphi + w0 * sinphi) / b; - uvw[1] = mu * v0 - a * b * cosphi; - uvw[2] = mu * w0 + a * (v0 * w0 * cosphi - u0 * sinphi) / b; + b = std::sqrt(1. - u.y*u.y); + return {mu*u.x + a*(u.x*u.y*cosphi + u.z*sinphi) / b, + mu*u.y - a*b*cosphi, + mu*u.z + a*(u.y*u.z*cosphi - u.x*sinphi) / b}; } } -Direction rotate_angle(Direction u, double mu, double* phi) -{ - double uvw[] {u.x, u.y, u.z}; - rotate_angle_c(uvw, mu, phi); - return {uvw[0], uvw[1], uvw[2]}; -} - - double maxwell_spectrum(double T) { // Set the random numbers double r1 = prn(); diff --git a/src/mesh.cpp b/src/mesh.cpp index dceddd9b8..97ec81ba2 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -385,9 +385,9 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, // just a bit for the purposes of determining if there was an intersection // in case the mesh surfaces coincide with lattice/geometric surfaces which // might produce finite-precision errors. - Position last_r {p->last_xyz}; - Position r {p->coord[0].xyz}; - Direction u {p->coord[0].uvw}; + Position last_r {p->r_last_}; + Position r {p->r()}; + Direction u {p->u()}; Position r0 = last_r + TINY_BIT*u; Position r1 = r - TINY_BIT*u; @@ -542,9 +542,9 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins // Determine if the track intersects the tally mesh. // Copy the starting and ending coordinates of the particle. - Position r0 {p->last_xyz_current}; - Position r1 {p->coord[0].xyz}; - Direction u {p->coord[0].uvw}; + Position r0 {p->r_last_current_}; + Position r1 {p->r()}; + Direction u {p->u()}; // Determine indices for starting and ending location. int n = n_dimension_; @@ -695,7 +695,7 @@ void RegularMesh::to_hdf5(hid_t group) const close_group(mesh_group); } -xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, +xt::xarray RegularMesh::count_sites(int64_t n, const Particle::Bank* bank, int n_energy, const double* energies, bool* outside) const { // Determine shape of array for counts @@ -713,7 +713,7 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, for (int64_t i = 0; i < n; ++i) { // determine scoring bin for entropy mesh - int mesh_bin = get_bin({bank[i].xyz}); + int mesh_bin = get_bin(bank[i].r); // if outside mesh, skip particle if (mesh_bin < 0) { @@ -771,7 +771,7 @@ openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end) { if (index_start) *index_start = model::meshes.size(); for (int i = 0; i < n; ++i) { - model::meshes.emplace_back(new RegularMesh{}); + model::meshes.push_back(std::make_unique()); } if (index_end) *index_end = model::meshes.size() - 1; @@ -909,7 +909,7 @@ void read_meshes(pugi::xml_node root) { for (auto node : root.children("mesh")) { // Read mesh and add to vector - model::meshes.emplace_back(new RegularMesh{node}); + model::meshes.push_back(std::make_unique(node)); // Map ID to position in vector model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1; diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 5b37d8ad5..d345cdb22 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -616,7 +616,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt) //============================================================================== void -Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3], +Mgxs::calculate_xs(int gin, double sqrtkT, Direction u, double& total_xs, double& abs_xs, double& nu_fiss_xs) { // Set our indices @@ -626,7 +626,7 @@ Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3], int tid = 0; #endif set_temperature_index(sqrtkT); - set_angle_index(uvw); + set_angle_index(u); XsData* xs_t = &xs[cache[tid].t]; total_xs = xs_t->total(cache[tid].a, gin); abs_xs = xs_t->absorption(cache[tid].a, gin); @@ -668,7 +668,7 @@ Mgxs::set_temperature_index(double sqrtkT) //============================================================================== void -Mgxs::set_angle_index(const double uvw[3]) +Mgxs::set_angle_index(Direction u) { // See if we need to find the new index #ifdef _OPENMP @@ -677,11 +677,11 @@ Mgxs::set_angle_index(const double uvw[3]) int tid = 0; #endif if (!is_isotropic && - ((uvw[0] != cache[tid].u) || (uvw[1] != cache[tid].v) || - (uvw[2] != cache[tid].w))) { - // convert uvw to polar and azimuthal angles - double my_pol = std::acos(uvw[2]); - double my_azi = std::atan2(uvw[1], uvw[0]); + ((u.x != cache[tid].u) || (u.y != cache[tid].v) || + (u.z != cache[tid].w))) { + // convert direction to polar and azimuthal angles + double my_pol = std::acos(u.z); + double my_azi = std::atan2(u.y, u.x); // Find the location, assuming equal-bin angles double delta_angle = PI / n_pol; @@ -692,9 +692,9 @@ Mgxs::set_angle_index(const double uvw[3]) cache[tid].a = n_azi * p + a; // store this direction as the last one used - cache[tid].u = uvw[0]; - cache[tid].v = uvw[1]; - cache[tid].w = uvw[2]; + cache[tid].u = u.x; + cache[tid].v = u.y; + cache[tid].w = u.z; } } diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index d2913a17a..c30850e8f 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -136,7 +136,7 @@ void create_macro_xs() for (int i = 0; i < model::materials.size(); ++i) { if (kTs[i].size() > 0) { // Convert atom_densities to a vector - Material* mat = model::materials[i]; + auto& mat {model::materials[i]}; std::vector atom_densities(mat->atom_density_.begin(), mat->atom_density_.end()); @@ -235,7 +235,7 @@ void read_mg_cross_sections_header() } // Get the minimum and maximum energies - int neutron = static_cast(ParticleType::neutron); + int neutron = static_cast(Particle::Type::neutron); data::energy_min[neutron] = data::energy_bins.back(); data::energy_max[neutron] = data::energy_bins.front(); @@ -248,10 +248,10 @@ void read_mg_cross_sections_header() //============================================================================== void -calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], +calculate_xs_c(int i_mat, int gin, double sqrtkT, Direction u, double& total_xs, double& abs_xs, double& nu_fiss_xs) { - data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs, + data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, u, total_xs, abs_xs, nu_fiss_xs); } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index a279abf2b..236833ace 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -287,7 +287,7 @@ void Nuclide::create_derived() auto xs = xt::adapt(rx->xs_[t].value); for (const auto& p : rx->products_) { - if (p.particle_ == ParticleType::photon) { + if (p.particle_ == Particle::Type::photon) { auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD); for (int k = 0; k < n; ++k) { double E = grid_[t].energy[k+j]; @@ -391,7 +391,7 @@ void Nuclide::create_derived() void Nuclide::init_grid() { - int neutron = static_cast(ParticleType::neutron); + int neutron = static_cast(Particle::Type::neutron); double E_min = data::energy_min[neutron]; double E_max = data::energy_max[neutron]; int M = settings::n_log_bins; @@ -440,7 +440,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const for (int i = 1; i < rx->products_.size(); ++i) { // Skip any non-neutron products const auto& product = rx->products_[i]; - if (product.particle_ != ParticleType::neutron) continue; + if (product.particle_ != Particle::Type::neutron) continue; // Evaluate yield if (product.emission_mode_ == EmissionMode::delayed) { diff --git a/src/output.cpp b/src/output.cpp index 095a17526..92ad14e9d 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -142,64 +142,64 @@ std::string time_stamp() extern "C" void print_particle(Particle* p) { // Display particle type and ID. - switch (p->type) { - case static_cast(ParticleType::neutron): + switch (p->type_) { + case Particle::Type::neutron: std::cout << "Neutron "; break; - case static_cast(ParticleType::photon): + case Particle::Type::photon: std::cout << "Photon "; break; - case static_cast(ParticleType::electron): + case Particle::Type::electron: std::cout << "Electron "; break; - case static_cast(ParticleType::positron): + case Particle::Type::positron: std::cout << "Positron "; break; default: std::cout << "Unknown Particle "; } - std::cout << p->id << "\n"; + std::cout << p->id_ << "\n"; // Display particle geometry hierarchy. - for (auto i = 0; i < p->n_coord; i++) { + for (auto i = 0; i < p->n_coord_; i++) { std::cout << " Level " << i << "\n"; - if (p->coord[i].cell != C_NONE) { - const Cell& c {*model::cells[p->coord[i].cell]}; + if (p->coord_[i].cell != C_NONE) { + const Cell& c {*model::cells[p->coord_[i].cell]}; std::cout << " Cell = " << c.id_ << "\n"; } - if (p->coord[i].universe != C_NONE) { - const Universe& u {*model::universes[p->coord[i].universe]}; + if (p->coord_[i].universe != C_NONE) { + const Universe& u {*model::universes[p->coord_[i].universe]}; std::cout << " Universe = " << u.id_ << "\n"; } - if (p->coord[i].lattice != F90_NONE) { - const Lattice& lat {*model::lattices[p->coord[i].lattice]}; + if (p->coord_[i].lattice != C_NONE) { + const Lattice& lat {*model::lattices[p->coord_[i].lattice]}; std::cout << " Lattice = " << lat.id_ << "\n"; - std::cout << " Lattice position = (" << p->coord[i].lattice_x - << "," << p->coord[i].lattice_y << "," - << p->coord[i].lattice_z << ")\n"; + std::cout << " Lattice position = (" << p->coord_[i].lattice_x + << "," << p->coord_[i].lattice_y << "," + << p->coord_[i].lattice_z << ")\n"; } - std::cout << " xyz = " << p->coord[i].xyz[0] << " " - << p->coord[i].xyz[1] << " " << p->coord[i].xyz[2] << "\n"; - std::cout << " uvw = " << p->coord[i].uvw[0] << " " - << p->coord[i].uvw[1] << " " << p->coord[i].uvw[2] << "\n"; + std::cout << " r = (" << p->coord_[i].r.x << ", " + << p->coord_[i].r.y << ", " << p->coord_[i].r.z << ")\n"; + std::cout << " u = (" << p->coord_[i].u.x << ", " + << p->coord_[i].u.y << ", " << p->coord_[i].u.z << ")\n"; } // Display miscellaneous info. - if (p->surface != ERROR_INT) { - const Surface& surf {*model::surfaces[std::abs(p->surface)-1]}; - std::cout << " Surface = " << std::copysign(surf.id_, p->surface) << "\n"; + if (p->surface_ != ERROR_INT) { + const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; + std::cout << " Surface = " << std::copysign(surf.id_, p->surface_) << "\n"; } - std::cout << " Weight = " << p->wgt << "\n"; + std::cout << " Weight = " << p->wgt_ << "\n"; if (settings::run_CE) { - std::cout << " Energy = " << p->E << "\n"; + std::cout << " Energy = " << p->E_ << "\n"; } else { - std::cout << " Energy Group = " << p->g << "\n"; + std::cout << " Energy Group = " << p->g_ << "\n"; } - std::cout << " Delayed Group = " << p->delayed_group << "\n"; + std::cout << " Delayed Group = " << p->delayed_group_ << "\n"; std::cout << "\n"; } diff --git a/src/particle.cpp b/src/particle.cpp index ccce1502f..ca7fccfe7 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -37,7 +37,7 @@ LocalCoord::reset() { cell = C_NONE; universe = C_NONE; - lattice = 0; + lattice = C_NONE; lattice_x = 0; lattice_y = 0; rotated = false; @@ -47,90 +47,70 @@ LocalCoord::reset() // Particle implementation //============================================================================== -void -Particle::clear() -{ - // reset any coordinate levels - for (int i=0; i(ParticleType::neutron); - alive = true; - - // clear attributes - surface = 0; - cell_born = C_NONE; - material = C_NONE; - last_material = C_NONE; - last_sqrtkT = 0; - wgt = 1.0; - last_wgt = 1.0; - absorb_wgt = 0.0; - n_bank = 0; - wgt_bank = 0.0; - sqrtkT = -1.0; - n_collision = 0; - fission = false; - delayed_group = 0; - for (int i=0; ir(); + secondary_bank_[n].u = u; + secondary_bank_[n].E = settings::run_CE ? E : g_; + ++n_secondary_; } void Particle::from_source(const Bank* src) { - // set defaults - initialize(); + // reset some attributes + this->clear(); + alive_ = true; + surface_ = 0; + cell_born_ = C_NONE; + material_ = C_NONE; + n_collision_ = 0; + fission_ = false; // copy attributes from source bank site - type = src->particle; - wgt = src->wgt; - last_wgt = src->wgt; - std::copy(src->xyz, src->xyz + 3, coord[0].xyz); - std::copy(src->uvw, src->uvw + 3, coord[0].uvw); - std::copy(src->xyz, src->xyz + 3, last_xyz_current); - std::copy(src->xyz, src->xyz + 3, last_xyz); - std::copy(src->uvw, src->uvw + 3, last_uvw); + type_ = src->particle; + wgt_ = src->wgt; + wgt_last_ = src->wgt; + this->r() = src->r; + this->u() = src->u; + r_last_current_ = src->r; + r_last_ = src->r; + u_last_ = src->u; if (settings::run_CE) { - E = src->E; - g = 0; + E_ = src->E; + g_ = 0; } else { - g = static_cast(src->E); - last_g = static_cast(src->E); - E = data::energy_bin_avg[g - 1]; + g_ = static_cast(src->E); + g_last_ = static_cast(src->E); + E_ = data::energy_bin_avg[g_ - 1]; } - last_E = E; + E_last_ = E_; } void @@ -138,7 +118,7 @@ Particle::transport() { // Display message if high verbosity or trace is on if (settings::verbosity >= 9 || simulation::trace) { - write_message("Simulating Particle " + std::to_string(id)); + write_message("Simulating Particle " + std::to_string(id_)); } // Initialize number of events to zero @@ -146,7 +126,7 @@ Particle::transport() // Add paricle's starting weight to count for normalizing tallies later #pragma omp atomic - simulation::total_weight += wgt; + simulation::total_weight += wgt_; // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { @@ -156,62 +136,62 @@ Particle::transport() } // Prepare to write out particle track. - if (write_track) add_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 == static_cast(ParticleType::neutron)) { + if (type_ == Particle::Type::neutron) { prn_set_stream(STREAM_TRACKING); } else { prn_set_stream(STREAM_PHOTON); } // Store pre-collision particle properties - last_wgt = wgt; - last_E = E; - std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw); - std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz); + wgt_last_ = wgt_; + E_last_ = E_; + u_last_ = this->u(); + r_last_ = this->r(); // 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 (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)); + + std::to_string(id_)); return; } // set birth cell attribute - if (cell_born == C_NONE) cell_born = coord[n_coord - 1].cell; + if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; } // Write particle track. - if (write_track) write_particle_track(*this); + 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 (material_ != MATERIAL_VOID) { if (settings::run_CE) { - if (material != last_material || sqrtkT != last_sqrtkT) { + 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); + model::materials[material_]->calculate_xs(*this); } } else { // Get the MG data - calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw, + calculate_xs_c(material_, g_, sqrtkT_, this->u_local(), simulation::material_xs.total, simulation::material_xs.absorption, simulation::material_xs.nu_fission); // Finally, update the particle group while we have already checked // for if multi-group - last_g = g; + g_last_ = g_; } } else { simulation::material_xs.total = 0.0; @@ -230,8 +210,8 @@ Particle::transport() // Sample a distance to collision double d_collision; - if (type == static_cast(ParticleType::electron) || - type == static_cast(ParticleType::positron)) { + if (type_ == Particle::Type::electron || + type_ == Particle::Type::positron) { d_collision = 0.0; } else if (simulation::material_xs.total == 0.0) { d_collision = INFINITY; @@ -243,11 +223,8 @@ Particle::transport() double distance = std::min(d_boundary, d_collision); // Advance particle - for (int j = 0; j < n_coord; ++j) { - // TODO: use Position - coord[j].xyz[0] += distance * coord[j].uvw[0]; - coord[j].xyz[1] += distance * coord[j].uvw[1]; - coord[j].xyz[2] += distance * coord[j].uvw[2]; + for (int j = 0; j < n_coord_; ++j) { + coord_[j].r += distance * coord_[j].u; } // Score track-length tallies @@ -257,8 +234,8 @@ Particle::transport() // Score track-length estimate of k-eff if (settings::run_mode == RUN_MODE_EIGENVALUE && - type == static_cast(ParticleType::neutron)) { - global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission; + type_ == Particle::Type::neutron) { + global_tally_tracklength += wgt_ * distance * simulation::material_xs.nu_fission; } // Score flux derivative accumulators for differential tallies. @@ -270,25 +247,25 @@ Particle::transport() // ==================================================================== // PARTICLE CROSSES SURFACE - if (next_level > 0) n_coord = next_level; + if (next_level > 0) n_coord_ = next_level; // Saving previous cell data - for (int j = 0; j < n_coord; ++j) { - last_cell[j] = coord[j].cell; + for (int j = 0; j < n_coord_; ++j) { + cell_last_[j] = coord_[j].cell; } - last_n_coord = n_coord; + n_coord_last_ = n_coord_; if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || lattice_translation[2] != 0) { // Particle crosses lattice boundary - surface = ERROR_INT; + surface_ = ERROR_INT; cross_lattice(this, lattice_translation); - event = EVENT_LATTICE; + event_ = EVENT_LATTICE; } else { // Particle crosses surface - surface = surface_crossed; + surface_ = surface_crossed; this->cross_surface(); - event = EVENT_SURFACE; + event_ = EVENT_SURFACE; } // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { @@ -300,8 +277,8 @@ Particle::transport() // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && - type == static_cast(ParticleType::neutron)) { - global_tally_collision += wgt * simulation::material_xs.nu_fission + type_ == Particle::Type::neutron) { + global_tally_collision += wgt_ * simulation::material_xs.nu_fission / simulation::material_xs.total; } @@ -313,7 +290,7 @@ Particle::transport() score_surface_tally(this, model::active_meshsurf_tallies); // Clear surface component - surface = ERROR_INT; + surface_ = ERROR_INT; if (settings::run_CE) { collision(this); @@ -334,33 +311,33 @@ Particle::transport() } // Reset banked weight during collision - n_bank = 0; - wgt_bank = 0.0; - for (int& v : n_delayed_bank) v = 0; + n_bank_ = 0; + wgt_bank_ = 0.0; + for (int& v : n_delayed_bank_) v = 0; // Reset fission logical - fission = false; + fission_ = false; // Save coordinates for tallying purposes - std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current); + r_last_current_ = this->r(); // Set last material to none since cross sections will need to be // re-evaluated - last_material = C_NONE; + material_last_ = C_NONE; - // Set all uvws to base level -- right now, after a collision, only the - // base level uvws are changed - for (int j = 0; j < n_coord - 1; ++j) { - if (coord[j + 1].rotated) { + // 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_}; - Direction u {coord[j].uvw}; - coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z; - coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z; - coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z; + const auto& m {model::cells[coord_[j].cell]->rotation_}; + const auto& u {coord_[j].u}; + coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z; + coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z; + coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z; } else { // Otherwise, copy this level's direction - std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw); + coord_[j+1].u = coord_[j].u; } } @@ -371,27 +348,27 @@ Particle::transport() // If particle has too many events, display warning and kill it ++n_event; if (n_event == MAX_EVENTS) { - warning("Particle " + std::to_string(id) + + warning("Particle " + std::to_string(id_) + " underwent maximum number of events."); - alive = false; + alive_ = false; } // Check for secondary particles if this particle is dead - if (!alive) { + if (!alive_) { // If no secondary particles, break out of event loop - if (n_secondary == 0) break; + if (n_secondary_ == 0) break; - this->from_source(&secondary_bank[n_secondary - 1]); - --n_secondary; + this->from_source(&secondary_bank_[n_secondary_ - 1]); + --n_secondary_; n_event = 0; // Enter new particle in particle track file - if (write_track) add_particle_track(); + if (write_track_) add_particle_track(); } } // Finish particle track output. - if (write_track) { + if (write_track_) { write_particle_track(*this); finalize_particle_track(*this); } @@ -400,9 +377,9 @@ Particle::transport() void Particle::cross_surface() { - int i_surface = std::abs(surface); + int i_surface = std::abs(surface_); // TODO: off-by-one - const auto& surf {model::surfaces[i_surface - 1]}; + const auto& surf {model::surfaces[i_surface - 1].get()}; if (settings::verbosity >= 10 || simulation::trace) { write_message(" Crossing surface " + std::to_string(surf->id_)); } @@ -412,7 +389,7 @@ Particle::cross_surface() // PARTICLE LEAKS OUT OF PROBLEM // Kill particle - alive = false; + alive_ = false; // Score any surface current tallies -- note that the particle is moved // forward slightly so that if the mesh boundary is on the surface, it is @@ -422,15 +399,12 @@ Particle::cross_surface() // TODO: Find a better solution to score surface currents than // physically moving the particle forward slightly - // TODO: Use Position - coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + this->r() += TINY_BIT * this->u(); score_surface_tally(this, model::active_meshsurf_tallies); } // Score to global leakage tally - global_tally_leakage += wgt; + global_tally_leakage += wgt_; // Display message if (settings::verbosity >= 10 || simulation::trace) { @@ -443,8 +417,8 @@ Particle::cross_surface() // PARTICLE REFLECTS FROM SURFACE // Do not handle reflective boundary conditions on lower universes - if (n_coord != 1) { - this->mark_as_lost("Cannot reflect particle " + std::to_string(id) + + if (n_coord_ != 1) { + this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) + " off surface in a lower universe."); return; } @@ -462,32 +436,27 @@ Particle::cross_surface() if (!model::active_meshsurf_tallies.empty()) { - Position r {coord[0].xyz}; - coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + Position r {this->r()}; + this->r() -= TINY_BIT * this->u(); score_surface_tally(this, model::active_meshsurf_tallies); - std::copy(&r.x, &r.x + 3, coord[0].xyz); + this->r() = r; } // Reflect particle off surface - Direction u = surf->reflect(coord[0].xyz, coord[0].uvw); + Direction u = surf->reflect(this->r(), this->u()); // Make sure new particle direction is normalized - double norm = u.norm(); - coord[0].uvw[0] = u.x/norm; - coord[0].uvw[1] = u.y/norm; - coord[0].uvw[2] = u.z/norm; + this->u() = u / u.norm(); // Reassign particle's cell and surface - coord[0].cell = last_cell[last_n_coord - 1]; - surface = -surface; + coord_[0].cell = cell_last_[n_coord_last_ - 1]; + surface_ = -surface_; // If a reflective surface is coincident with a lattice or universe // boundary, it is necessary to redetermine the particle's coordinates in // the lower universes. - n_coord = 1; + n_coord_ = 1; if (!find_cell(this, true)) { this->mark_as_lost("Couldn't find particle after reflecting from surface " + std::to_string(surf->id_) + "."); @@ -495,9 +464,7 @@ Particle::cross_surface() } // Set previous coordinate going slightly past surface crossing - last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT*coord[0].uvw[0]; - last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT*coord[0].uvw[1]; - last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT*coord[0].uvw[2]; + r_last_current_ = this->r() + TINY_BIT*this->u(); // Diagnostic message if (settings::verbosity >= 10 || simulation::trace) { @@ -510,8 +477,8 @@ Particle::cross_surface() // PERIODIC BOUNDARY // Do not handle periodic boundary conditions on lower universes - if (n_coord != 1) { - this->mark_as_lost("Cannot transfer particle " + std::to_string(id) + + if (n_coord_ != 1) { + this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) + " across surface in a lower universe. Boundary conditions must be " "applied to root universe."); return; @@ -521,34 +488,28 @@ Particle::cross_surface() // particle to change -- artificially move the particle slightly back in // case the surface crossing is coincident with a mesh boundary if (!model::active_meshsurf_tallies.empty()) { - Position r {coord[0].xyz}; - coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + Position r {this->r()}; + this->r() -= TINY_BIT * this->u(); score_surface_tally(this, model::active_meshsurf_tallies); - std::copy(&r.x, &r.x + 3, coord[0].xyz); + this->r() = r; } // Get a pointer to the partner periodic surface auto surf_p = dynamic_cast(surf); auto other = dynamic_cast( - model::surfaces[surf_p->i_periodic_]); + model::surfaces[surf_p->i_periodic_].get()); // Adjust the particle's location and direction. - Position r {coord[0].xyz}; - Direction u {coord[0].uvw}; - bool rotational = other->periodic_translate(surf_p, r, u); - std::copy(&r.x, &r.x + 3, coord[0].xyz); - std::copy(&u.x, &u.x + 3, coord[0].uvw); + bool rotational = other->periodic_translate(surf_p, this->r(), this->u()); // Reassign particle's surface // TODO: off-by-one - surface = rotational ? + surface_ = rotational ? surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface); + std::copysign(surf_p->i_periodic_ + 1, surface_); // Figure out what cell particle is in now - n_coord = 1; + n_coord_ = 1; if (!find_cell(this, true)) { this->mark_as_lost("Couldn't find particle after hitting periodic " @@ -557,9 +518,7 @@ Particle::cross_surface() } // Set previous coordinate going slightly past surface crossing - last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT * coord[0].uvw[0]; - last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT * coord[0].uvw[1]; - last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT * coord[0].uvw[2]; + r_last_current_ = this->r() + TINY_BIT*this->u(); // Diagnostic message if (settings::verbosity >= 10 || simulation::trace) { @@ -574,18 +533,18 @@ Particle::cross_surface() #ifdef DAGMC if (settings::dagmc) { - auto cellp = dynamic_cast(model::cells[last_cell[0]]); + auto cellp = dynamic_cast(model::cells[cell_last_[0]].get()); // TODO: off-by-one - auto surfp = dynamic_cast(model::surfaces[std::abs(surface) - 1]); + auto surfp = dynamic_cast(model::surfaces[std::abs(surface_) - 1].get()); int32_t i_cell = next_cell(cellp, surfp) - 1; // save material and temp - last_material = material; - last_sqrtkT = sqrtkT; + material_last_ = material_; + sqrtkT_last_ = sqrtkT_; // set new cell value - coord[0].cell = i_cell; - cell_instance = 0; - material = model::cells[i_cell]->material_[0]; - sqrtkT = model::cells[i_cell]->sqrtkT_[0]; + coord_[0].cell = i_cell; + cell_instance_ = 0; + material_ = model::cells[i_cell]->material_[0]; + sqrtkT_ = model::cells[i_cell]->sqrtkT_[0]; return; } #endif @@ -596,8 +555,8 @@ Particle::cross_surface() // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS // Remove lower coordinate levels and assignment of surface - surface = ERROR_INT; - n_coord = 1; + surface_ = ERROR_INT; + n_coord_ = 1; bool found = find_cell(this, false); if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) { @@ -606,16 +565,14 @@ Particle::cross_surface() // the particle is really traveling tangent to a surface, if we move it // forward a tiny bit it should fix the problem. - n_coord = 1; - coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + n_coord_ = 1; + this->r() += TINY_BIT * this->u(); // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. if (!find_cell(this, false)) { - this->mark_as_lost("After particle " + std::to_string(id) + + this->mark_as_lost("After particle " + std::to_string(id_) + " crossed surface " + std::to_string(surf->id_) + " it could not be located in any cell and it did not leak."); return; @@ -631,7 +588,7 @@ Particle::mark_as_lost(const char* message) write_restart(); // Increment number of lost particles - alive = false; + alive_ = false; #pragma omp atomic simulation::n_lost_particles += 1; @@ -655,7 +612,7 @@ Particle::write_restart() const // Set up file name std::stringstream filename; filename << settings::path_output << "particle_" << simulation::current_batch - << '_' << id << ".h5"; + << '_' << id_ << ".h5"; #pragma omp critical (WriteParticleRestart) { @@ -686,15 +643,14 @@ Particle::write_restart() const write_dataset(file_id, "run_mode", "particle restart"); break; } - write_dataset(file_id, "id", id); - write_dataset(file_id, "type", type); + write_dataset(file_id, "id", id_); + write_dataset(file_id, "type", static_cast(type_)); int64_t i = simulation::current_work; write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt); write_dataset(file_id, "energy", simulation::source_bank[i-1].E); - hsize_t dims[] {3}; - write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false); - write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false); + write_dataset(file_id, "xyz", simulation::source_bank[i-1].r); + write_dataset(file_id, "uvw", simulation::source_bank[i-1].u); // Close file file_close(file_id); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 19dd50d77..0cface319 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -39,29 +39,28 @@ void read_particle_restart(Particle& p, int& previous_run_mode) } else if (mode == "fixed source") { previous_run_mode = RUN_MODE_FIXEDSOURCE; } - read_dataset(file_id, "id", p.id); - read_dataset(file_id, "type", p.type); - read_dataset(file_id, "weight", p.wgt); - read_dataset(file_id, "energy", p.E); - std::array x; - read_dataset(file_id, "xyz", x); - std::copy(x.data(), x.data() + 3, p.coord[0].xyz); - read_dataset(file_id, "uvw", x); - std::copy(x.data(), x.data() + 3, p.coord[0].uvw); + read_dataset(file_id, "id", p.id_); + int type; + read_dataset(file_id, "type", type); + p.type_ = static_cast(type); + read_dataset(file_id, "weight", p.wgt_); + read_dataset(file_id, "energy", p.E_); + read_dataset(file_id, "xyz", p.r()); + read_dataset(file_id, "uvw", p.u()); // Set energy group and average energy in multi-group mode if (!settings::run_CE) { - p.g = p.E; - p.E = data::energy_bin_avg[p.g - 1]; + p.g_ = p.E_; + p.E_ = data::energy_bin_avg[p.g_ - 1]; } // Set particle last attributes - p.last_wgt = p.wgt; - std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current); - std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz); - std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw); - p.last_E = p.E; - p.last_g = p.g; + p.wgt_last_ = p.wgt_; + p.r_last_current_ = p.r(); + p.r_last_ = p.r(); + p.u_last_ = p.u(); + p.E_last_ = p.E_; + p.g_last_ = p.g_; // Close hdf5 file file_close(file_id); @@ -81,7 +80,6 @@ void run_particle_restart() // Initialize the particle to be tracked Particle p; - p.initialize(); // Read in the restart information int previous_run_mode; @@ -94,10 +92,10 @@ void run_particle_restart() int64_t particle_seed; switch (previous_run_mode) { case RUN_MODE_EIGENVALUE: - particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id; + particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id_; break; case RUN_MODE_FIXEDSOURCE: - particle_seed = p.id; + particle_seed = p.id_; break; } set_particle_seed(particle_seed); diff --git a/src/photon.cpp b/src/photon.cpp index 89fef6c43..2dd0928ba 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -224,7 +224,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) } // Truncate the bremsstrahlung data at the cutoff energy - int photon = static_cast(ParticleType::photon); + int photon = static_cast(Particle::Type::photon); const auto& E {electron_energy}; double cutoff = settings::energy_cutoff[photon]; if (cutoff > E(0)) { @@ -653,13 +653,12 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl if (shell.n_transitions == 0) { double mu = 2.0*prn() - 1.0; double phi = 2.0*PI*prn(); - std::array uvw; - uvw[0] = mu; - uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); - uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); double E = shell.binding_energy; - int photon = static_cast(ParticleType::photon); - p.create_secondary(uvw.data(), E, photon, true); + p.create_secondary(u, E, Particle::Type::photon); return; } @@ -679,10 +678,10 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Sample angle isotropically double mu = 2.0*prn() - 1.0; double phi = 2.0*PI*prn(); - std::array uvw; - uvw[0] = mu; - uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); - uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Get the transition energy double E = shell.transition_energy(i_transition); @@ -691,8 +690,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-radiative transition -- Auger/Coster-Kronig effect // Create auger electron - int electron = static_cast(ParticleType::electron); - p.create_secondary(uvw.data(), E, electron, true); + p.create_secondary(u, E, Particle::Type::electron); // Fill hole left by emitted auger electron int i_hole = shell_map_.at(secondary); @@ -702,8 +700,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Radiative transition -- get X-ray energy // Create fluorescent photon - int photon = static_cast(ParticleType::photon); - p.create_secondary(uvw.data(), E, photon, true); + p.create_secondary(u, E, Particle::Type::photon); } // Fill hole created by electron transitioning to the photoelectron hole diff --git a/src/physics.cpp b/src/physics.cpp index 5d6ca2dee..d898ed788 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -33,40 +33,41 @@ namespace openmc { void collision(Particle* p) { // Add to collision counter for particle - ++(p->n_collision); + ++(p->n_collision_); // Sample reaction for the material the particle is in - switch (static_cast(p->type)) { - case ParticleType::neutron: + switch (static_cast(p->type_)) { + case Particle::Type::neutron: sample_neutron_reaction(p); break; - case ParticleType::photon: + case Particle::Type::photon: sample_photon_reaction(p); break; - case ParticleType::electron: + case Particle::Type::electron: sample_electron_reaction(p); break; - case ParticleType::positron: + case Particle::Type::positron: sample_positron_reaction(p); break; } // Kill particle if energy falls below cutoff - if (p->E < settings::energy_cutoff[p->type]) { - p->alive = false; - p->wgt = 0.0; - p->last_wgt = 0.0; + int type = static_cast(p->type_); + if (p->E_ < settings::energy_cutoff[type]) { + p->alive_ = false; + p->wgt_ = 0.0; + p->wgt_last_ = 0.0; } // Display information about collision if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; - if (static_cast(p->type) == ParticleType::neutron) { - msg << " " << reaction_name(p->event_MT) << " with " << - data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV."; + if (p->type_ == Particle::Type::neutron) { + msg << " " << reaction_name(p->event_mt_) << " with " << + data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else { - msg << " " << reaction_name(p->event_MT) << " with " << - data::elements[p->event_nuclide].name_ << ". Energy = " << p->E << " eV."; + msg << " " << reaction_name(p->event_mt_) << " with " << + data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV."; } write_message(msg, 1); } @@ -78,7 +79,7 @@ void sample_neutron_reaction(Particle* p) int i_nuclide = sample_nuclide(p); // Save which nuclide particle had collision with - p->event_nuclide = i_nuclide; + p->event_nuclide_ = i_nuclide; // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not @@ -88,14 +89,14 @@ void sample_neutron_reaction(Particle* p) const auto& nuc {data::nuclides[i_nuclide]}; if (nuc->fissionable_) { - Reaction* rx = sample_fission(i_nuclide, p->E); + Reaction* rx = sample_fission(i_nuclide, p->E_); if (settings::run_mode == RUN_MODE_EIGENVALUE) { create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(), &simulation::n_bank, simulation::fission_bank.size()); } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx, p->secondary_bank, - &p->n_secondary, MAX_SECONDARY); + create_fission_sites(p, i_nuclide, rx, p->secondary_bank_, + &p->n_secondary_, MAX_SECONDARY); } } @@ -112,16 +113,16 @@ void sample_neutron_reaction(Particle* p) if (simulation::micro_xs[i_nuclide].absorption > 0.0) { absorption(p, i_nuclide); } else { - p->absorb_wgt = 0.0; + p->wgt_absorb_ = 0.0; } - if (!p->alive) return; + if (!p->alive_) return; // Sample a scattering reaction and determine the secondary energy of the // exiting neutron scatter(p, i_nuclide); // Advance URR seed stream 'N' times after energy changes - if (p->E != p->last_E) { + if (p->E_ != p->E_last_) { prn_set_stream(STREAM_URR_PTABLE); advance_prn_seed(data::nuclides.size()); prn_set_stream(STREAM_TRACKING); @@ -130,13 +131,13 @@ void sample_neutron_reaction(Particle* p) // Play russian roulette if survival biasing is turned on if (settings::survival_biasing) { russian_roulette(p); - if (!p->alive) return; + if (!p->alive_) return; } } void -create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array, - int64_t* size_bank, int64_t bank_capacity) +create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, + Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity) { // TODO: Heat generation from fission @@ -145,7 +146,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt / simulation::keff * weight * simulation::micro_xs[ + double nu_t = p->wgt_ / simulation::keff * weight * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].total; // Sample the number of neutrons produced @@ -177,28 +178,26 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; - p->fission = true; + p->fission_ = true; for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) { // Bank source neutrons by copying the particle data - bank_array[i].xyz[0] = p->coord[0].xyz[0]; - bank_array[i].xyz[1] = p->coord[0].xyz[1]; - bank_array[i].xyz[2] = p->coord[0].xyz[2]; + bank_array[i].r = p->r(); // Set that the bank particle is a neutron - bank_array[i].particle = static_cast(ParticleType::neutron); + bank_array[i].particle = Particle::Type::neutron; // Set the weight of the fission bank site bank_array[i].wgt = 1. / weight; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p->E, &bank_array[i]); + sample_fission_neutron(i_nuclide, rx, p->E_, &bank_array[i]); // Set the delayed group on the particle as well - p->delayed_group = bank_array[i].delayed_group; + p->delayed_group_ = bank_array[i].delayed_group; // Increment the number of neutrons born delayed - if (p->delayed_group > 0) { - nu_d[p->delayed_group-1]++; + if (p->delayed_group_ > 0) { + nu_d[p->delayed_group_-1]++; } } @@ -206,10 +205,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ *size_bank = std::min(*size_bank + nu, bank_capacity); // Store the total weight banked for analog fission tallies - p->n_bank = nu; - p->wgt_bank = nu / weight; + p->n_bank_ = nu; + p->wgt_bank_ = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p->n_delayed_bank[d] = nu_d[d]; + p->n_delayed_bank_[d] = nu_d[d]; } } @@ -218,21 +217,21 @@ void sample_photon_reaction(Particle* p) // Kill photon if below energy cutoff -- an extra check is made here because // photons with energy below the cutoff may have been produced by neutrons // reactions or atomic relaxation - int photon = static_cast(ParticleType::photon); - if (p->E < settings::energy_cutoff[photon]) { - p->E = 0.0; - p->alive = false; + int photon = static_cast(Particle::Type::photon); + if (p->E_ < settings::energy_cutoff[photon]) { + p->E_ = 0.0; + p->alive_ = false; return; } // Sample element within material int i_element = sample_element(p); - p->event_nuclide = i_element; + p->event_nuclide_ = i_element; const auto& micro {simulation::micro_photon_xs[i_element]}; const auto& element {data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent - double alpha = p->E/MASS_ELECTRON_EV; + double alpha = p->E_/MASS_ELECTRON_EV; // For tallying purposes, this routine might be called directly. In that // case, we need to sample a reaction via the cutoff variable @@ -243,8 +242,8 @@ void sample_photon_reaction(Particle* p) prob += micro.coherent; if (prob > cutoff) { double mu = element.rayleigh_scatter(alpha); - rotate_angle_c(p->coord[0].uvw, mu, nullptr); - p->event_MT = COHERENT; + p->u() = rotate_angle(p->u(), mu, nullptr); + p->event_mt_ = COHERENT; return; } @@ -269,11 +268,8 @@ void sample_photon_reaction(Particle* p) double mu_electron = (alpha - alpha_out*mu) / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); double phi = 2.0*PI*prn(); - double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); - rotate_angle_c(uvw, mu_electron, &phi); - int electron = static_cast(ParticleType::electron); - p->create_secondary(uvw, E_electron, electron, true); + Direction u = rotate_angle(p->u(), mu_electron, &phi); + p->create_secondary(u, E_electron, Particle::Type::electron); // TODO: Compton subshell data does not match atomic relaxation data // Allow electrons to fill orbital and produce auger electrons @@ -284,9 +280,9 @@ void sample_photon_reaction(Particle* p) } phi += PI; - p->E = alpha_out*MASS_ELECTRON_EV; - rotate_angle_c(p->coord[0].uvw, mu, &phi); - p->event_MT = INCOHERENT; + p->E_ = alpha_out*MASS_ELECTRON_EV; + p->u() = rotate_angle(p->u(), mu, &phi); + p->event_mt_ = INCOHERENT; return; } @@ -309,7 +305,7 @@ void sample_photon_reaction(Particle* p) prob += xs; if (prob > cutoff) { - double E_electron = p->E - shell.binding_energy; + double E_electron = p->E_ - shell.binding_energy; // Sample mu using non-relativistic Sauter distribution. // See Eqns 3.19 and 3.20 in "Implementing a photon physics @@ -326,21 +322,20 @@ void sample_photon_reaction(Particle* p) } double phi = 2.0*PI*prn(); - std::array uvw; - uvw[0] = mu; - uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); - uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create secondary electron - int electron = static_cast(ParticleType::electron); - p->create_secondary(uvw.data(), E_electron, electron, true); + p->create_secondary(u, E_electron, Particle::Type::electron); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons element.atomic_relaxation(shell, *p); - p->event_MT = 533 + shell.index_subshell; - p->alive = false; - p->E = 0.0; + p->event_mt_ = 533 + shell.index_subshell; + p->alive_ = false; + p->E_ = 0.0; return; } } @@ -356,21 +351,16 @@ void sample_photon_reaction(Particle* p) &mu_electron, &mu_positron); // Create secondary electron - double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); - rotate_angle_c(uvw, mu_electron, nullptr); - int electron = static_cast(ParticleType::electron); - p->create_secondary(uvw, E_electron, electron, true); + Direction u = rotate_angle(p->u(), mu_electron, nullptr); + p->create_secondary(u, E_electron, Particle::Type::electron); // Create secondary positron - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); - rotate_angle_c(uvw, mu_positron, nullptr); - int positron = static_cast(ParticleType::positron); - p->create_secondary(uvw, E_positron, positron, true); + u = rotate_angle(p->u(), mu_positron, nullptr); + p->create_secondary(u, E_positron, Particle::Type::positron); - p->event_MT = PAIR_PROD; - p->alive = false; - p->E = 0.0; + p->event_mt_ = PAIR_PROD; + p->alive_ = false; + p->E_ = 0.0; } } @@ -383,8 +373,8 @@ void sample_electron_reaction(Particle* p) thick_target_bremsstrahlung(*p, &E_lost); } - p->E = 0.0; - p->alive = false; + p->E_ = 0.0; + p->alive_ = false; } void sample_positron_reaction(Particle* p) @@ -399,22 +389,17 @@ void sample_positron_reaction(Particle* p) // Sample angle isotropically double mu = 2.0*prn() - 1.0; double phi = 2.0*PI*prn(); - std::array uvw; - uvw[0] = mu; - uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi); - uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create annihilation photon pair traveling in opposite directions - int photon = static_cast(ParticleType::photon); - p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true); + p->create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon); + p->create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon); - uvw[0] = -uvw[0]; - uvw[1] = -uvw[1]; - uvw[2] = -uvw[2]; - p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true); - - p->E = 0.0; - p->alive = false; + p->E_ = 0.0; + p->alive_ = false; } int sample_nuclide(const Particle* p) @@ -423,7 +408,7 @@ int sample_nuclide(const Particle* p) double cutoff = prn() * simulation::material_xs.total; // Get pointers to nuclide/density arrays - const auto& mat {model::materials[p->material]}; + const auto& mat {model::materials[p->material_]}; int n = mat->nuclide_.size(); double prob = 0.0; @@ -448,7 +433,7 @@ int sample_element(Particle* p) double cutoff = prn() * simulation::material_xs.total; // Get pointers to elements, densities - const auto& mat {model::materials[p->material]}; + const auto& mat {model::materials[p->material_]}; int n = mat->nuclide_.size(); int i = 0; @@ -541,7 +526,7 @@ void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product) + f*(rx->xs_[i_temp].value[i_grid - threshold + 1])); for (int j = 0; j < rx->products_.size(); ++j) { - if (rx->products_[j].particle_ == ParticleType::photon) { + if (rx->products_[j].particle_ == Particle::Type::photon) { // add to cumulative probability prob += (*rx->products_[j].yield_)(E) * xs; @@ -557,16 +542,16 @@ void absorption(Particle* p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * simulation::micro_xs[i_nuclide].absorption / + p->wgt_absorb_ = p->wgt_ * simulation::micro_xs[i_nuclide].absorption / simulation::micro_xs[i_nuclide].total; // Adjust weight of particle by probability of absorption - p->wgt -= p->absorb_wgt; - p->last_wgt = p->wgt; + p->wgt_ -= p->wgt_absorb_; + p->wgt_last_ = p->wgt_; // Score implicit absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->absorb_wgt * simulation::micro_xs[ + global_tally_absorption += p->wgt_absorb_ * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; } } else { @@ -575,13 +560,13 @@ void absorption(Particle* p, int i_nuclide) prn() * simulation::micro_xs[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt * simulation::micro_xs[ + global_tally_absorption += p->wgt_ * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; } - p->alive = false; - p->event = EVENT_ABSORB; - p->event_MT = N_DISAPPEAR; + p->alive_ = false; + p->event_ = EVENT_ABSORB; + p->event_mt_ = N_DISAPPEAR; } } } @@ -589,7 +574,7 @@ void absorption(Particle* p, int i_nuclide) void scatter(Particle* p, int i_nuclide) { // copy incoming direction - Direction u_old {p->coord[0].uvw}; + Direction u_old {p->u()}; // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; @@ -614,13 +599,13 @@ void scatter(Particle* p, int i_nuclide) // NON-S(A,B) ELASTIC SCATTERING // Determine temperature - double kT = nuc->multipole_ ? p->sqrtkT*p->sqrtkT : nuc->kTs_[i_temp]; + double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp]; // Perform collision physics for elastic scattering - elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT, - &p->E, p->coord[0].uvw, &p->mu, &p->wgt); + elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, + p->E_, p->u(), p->mu_); - p->event_MT = ELASTIC; + p->event_mt_ = ELASTIC; sampled = true; } @@ -629,9 +614,9 @@ void scatter(Particle* p, int i_nuclide) // ======================================================================= // S(A,B) SCATTERING - sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu); + sab_scatter(i_nuclide, micro.index_sab, p->E_, p->u(), p->mu_); - p->event_MT = ELASTIC; + p->event_mt_ = ELASTIC; sampled = true; } @@ -663,53 +648,47 @@ void scatter(Particle* p, int i_nuclide) // Perform collision physics for inelastic scattering const auto& rx {nuc->reactions_[i]}; inelastic_scatter(nuc.get(), rx.get(), p); - p->event_MT = rx->mt_; + p->event_mt_ = rx->mt_; } // Set event component - p->event = EVENT_SCATTER; + p->event_ = EVENT_SCATTER; // Sample new outgoing angle for isotropic-in-lab scattering - const auto& mat {model::materials[p->material]}; + const auto& mat {model::materials[p->material_]}; if (!mat->p0_.empty()) { int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; if (mat->p0_[i_nuc_mat]) { // Sample isotropic-in-lab outgoing direction double mu = 2.0*prn() - 1.0; double phi = 2.0*PI*prn(); - Direction u_new; - u_new.x = mu; - u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); - u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); - - p->mu = u_old.dot(u_new); // Change direction of particle - p->coord[0].uvw[0] = u_new.x; - p->coord[0].uvw[1] = u_new.y; - p->coord[0].uvw[2] = u_new.z; + p->u().x = mu; + p->u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + p->u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi); + p->mu_ = u_old.dot(p->u()); } } } -void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, - double* uvw, double* mu_lab, double* wgt) +void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E, + Direction& u, double& mu_lab) { // get pointer to nuclide const auto& nuc {data::nuclides[i_nuclide]}; - double vel = std::sqrt(*E); + double vel = std::sqrt(E); double awr = nuc->awr_; // Neutron velocity in LAB - Direction u {uvw}; Direction v_n = vel*u; // Sample velocity of target nucleus Direction v_t {}; if (!simulation::micro_xs[i_nuclide].use_ptable) { - v_t = sample_target_velocity(nuc.get(), *E, u, v_n, - simulation::micro_xs[i_nuclide].elastic, kT, wgt); + v_t = sample_target_velocity(nuc.get(), E, u, v_n, + simulation::micro_xs[i_nuclide].elastic, kT); } // Velocity of center-of-mass @@ -724,10 +703,10 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, // Sample scattering angle, checking if it is an ncorrelated angle-energy // distribution double mu_cm; - auto& d = rx->products_[0].distribution_[0]; + auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); if (d_) { - mu_cm = d_->angle().sample(*E); + mu_cm = d_->angle().sample(E); } else { mu_cm = 2.0*prn() - 1.0; } @@ -743,26 +722,23 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, // Transform back to LAB frame v_n += v_cm; - *E = v_n.dot(v_n); - vel = std::sqrt(*E); + E = v_n.dot(v_n); + vel = std::sqrt(E); // compute cosine of scattering angle in LAB frame by taking dot product of // neutron's pre- and post-collision angle - *mu_lab = u.dot(v_n) / vel; + mu_lab = u.dot(v_n) / vel; // Set energy and direction of particle in LAB frame u = v_n / vel; - uvw[0] = u.x; - uvw[1] = u.y; - uvw[2] = u.z; // Because of floating-point roundoff, it may be possible for mu_lab to be // outside of the range [-1,1). In these cases, we just set mu_lab to exactly // -1 or 1 - if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab); + if (std::abs(mu_lab) > 1.0) mu_lab = std::copysign(1.0, mu_lab); } -void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu) +void sab_scatter(int i_nuclide, int i_sab, double& E, Direction& u, double& mu) { // Determine temperature index const auto& micro {simulation::micro_xs[i_nuclide]}; @@ -770,15 +746,15 @@ void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu) // Sample energy and angle double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, *E, &E_out, mu); + data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, E, &E_out, &mu); // Set energy to outgoing, change direction of particle - *E = E_out; - rotate_angle_c(uvw, *mu, nullptr); + E = E_out; + u = rotate_angle(u, mu, nullptr); } Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, double* wgt) + Direction v_neut, double xs_eff, double kT) { // check if nuclide is a resonant scatterer ResScatMethod sampling_method; @@ -977,7 +953,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT) return vt * rotate_angle(u, mu, nullptr); } -void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site) +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site) { // Sample cosine of angle -- fission neutrons are always emitted // isotropically. Sometimes in ACE data, fission reactions actually have @@ -987,9 +963,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank // Sample azimuthal angle uniformly in [0,2*pi) double phi = 2.0*PI*prn(); - site->uvw[0] = mu; - site->uvw[1] = std::sqrt(1.0 - mu*mu) * std::cos(phi); - site->uvw[2] = std::sqrt(1.0 - mu*mu) * std::sin(phi); + site->u.x = mu; + site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi); + site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi); // Determine total nu, delayed nu, and delayed neutron fraction const auto& nuc {data::nuclides[i_nuclide]}; @@ -1029,7 +1005,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank rx->products_[group].sample(E_in, site->E, mu); // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(ParticleType::neutron); + constexpr int neutron = static_cast(Particle::Type::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples @@ -1054,7 +1030,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank rx->products_[0].sample(E_in, site->E, mu); // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(ParticleType::neutron); + constexpr int neutron = static_cast(Particle::Type::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples @@ -1071,7 +1047,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) { // copy energy of neutron - double E_in = p->E; + double E_in = p->E_; // sample outgoing energy and scattering cosine double E; @@ -1098,30 +1074,29 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); // Set outgoing energy and scattering angle - p->E = E; - p->mu = mu; + p->E_ = E; + p->mu_ = mu; // change direction of particle - rotate_angle_c(p->coord[0].uvw, mu, nullptr); + p->u() = rotate_angle(p->u(), mu, nullptr); // evaluate yield double yield = (*rx->products_[0].yield_)(E_in); if (std::floor(yield) == yield) { // If yield is integral, create exactly that many secondary particles for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - int neutron = static_cast(ParticleType::neutron); - p->create_secondary(p->coord[0].uvw, p->E, neutron, true); + p->create_secondary(p->u(), p->E_, Particle::Type::neutron); } } else { // Otherwise, change weight of particle based on yield - p->wgt *= yield; + p->wgt_ *= yield; } } void sample_secondary_photons(Particle* p, int i_nuclide) { // Sample the number of photons produced - double y_t = p->wgt * simulation::micro_xs[i_nuclide].photon_prod / + double y_t = p->wgt_ * simulation::micro_xs[i_nuclide].photon_prod / simulation::micro_xs[i_nuclide].total; int y = static_cast(y_t); if (prn() <= y_t - y) ++y; @@ -1131,22 +1106,19 @@ void sample_secondary_photons(Particle* p, int i_nuclide) // Sample the reaction and product int i_rx; int i_product; - sample_photon_product(i_nuclide, p->E, &i_rx, &i_product); + sample_photon_product(i_nuclide, p->E_, &i_rx, &i_product); // Sample the outgoing energy and angle auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; double E; double mu; - rx->products_[i_product].sample(p->E, E, mu); + rx->products_[i_product].sample(p->E_, E, mu); // Sample the new direction - double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); - rotate_angle_c(uvw, mu, nullptr); + Direction u = rotate_angle(p->u(), mu, nullptr); // Create the secondary photon - int photon = static_cast(ParticleType::photon); - p->create_secondary(uvw, E, photon, true); + p->create_secondary(u, E, Particle::Type::photon); } } diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 02c77ae7d..42ddc2ecb 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -11,14 +11,14 @@ namespace openmc { void russian_roulette(Particle* p) { - if (p->wgt < settings::weight_cutoff) { - if (prn() < p->wgt / settings::weight_survive) { - p->wgt = settings::weight_survive; - p->last_wgt = p->wgt; + if (p->wgt_ < settings::weight_cutoff) { + if (prn() < p->wgt_ / settings::weight_survive) { + p->wgt_ = settings::weight_survive; + p->wgt_last_ = p->wgt_; } else { - p->wgt = 0.; - p->last_wgt = 0.; - p->alive = false; + p->wgt_ = 0.; + p->wgt_last_ = 0.; + p->alive_ = false; } } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 87dd37ef8..fa6365a3a 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -25,7 +25,7 @@ void collision_mg(Particle* p) { // Add to the collision counter for the particle - p->n_collision++; + p->n_collision_++; // Sample the reaction type sample_reaction(p); @@ -33,7 +33,7 @@ collision_mg(Particle* p) // Display information about collision if ((settings::verbosity >= 10) || (simulation::trace)) { std::stringstream msg; - msg << " Energy Group = " << p->g; + msg << " Energy Group = " << p->g_; write_message(msg, 1); } } @@ -46,14 +46,14 @@ sample_reaction(Particle* p) // change when sampling fission sites. The following block handles all // absorption (including fission) - if (model::materials[p->material]->fissionable_) { + if (model::materials[p->material_]->fissionable_) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { create_fission_sites( p, simulation::fission_bank.data(), &simulation::n_bank, simulation::fission_bank.size()); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { - create_fission_sites(p, p->secondary_bank, &(p->n_secondary), + create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_), MAX_SECONDARY); } } @@ -63,9 +63,9 @@ sample_reaction(Particle* p) if (simulation::material_xs.absorption > 0.) { absorption(p); } else { - p->absorb_wgt = 0.; + p->wgt_absorb_ = 0.; } - if (!p->alive) return; + if (!p->alive_) return; // Sample a scattering event to determine the energy of the exiting neutron scatter(p); @@ -73,7 +73,7 @@ sample_reaction(Particle* p) // Play Russian roulette if survival biasing is turned on if (settings::survival_biasing) { russian_roulette(p); - if (!p->alive) return; + if (!p->alive_) return; } } @@ -82,27 +82,27 @@ scatter(Particle* p) { // Adjust indices for Fortran to C++ indexing // TODO: Remove when no longer needed - int gin = p->last_g - 1; - int gout = p->g - 1; - int i_mat = p->material; - data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt); + int gin = p->g_last_ - 1; + int gout = p->g_ - 1; + int i_mat = p->material_; + data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_); // Adjust return value for fortran indexing // TODO: Remove when no longer needed - p->g = gout + 1; + p->g_ = gout + 1; // Rotate the angle - rotate_angle_c(p->coord[0].uvw, p->mu, nullptr); + p->u() = rotate_angle(p->u(), p->mu_, nullptr); // Update energy value for downstream compatability (in tallying) - p->E = data::energy_bin_avg[gout]; + p->E_ = data::energy_bin_avg[gout]; // Set event component - p->event = EVENT_SCATTER; + p->event_ = EVENT_SCATTER; } void -create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, +create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_array_size) { // TODO: Heat generation from fission @@ -112,7 +112,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt / simulation::keff * weight * + double nu_t = p->wgt_ / simulation::keff * weight * simulation::material_xs.nu_fission / simulation::material_xs.total; // Sample the number of neutrons produced @@ -148,16 +148,14 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; - p->fission = true; + p->fission_ = true; for (size_t i = static_cast(*size_bank); i < static_cast(std::min(*size_bank + nu, bank_array_size)); i++) { // Bank source neutrons by copying the particle data - bank_array[i].xyz[0] = p->coord[0].xyz[0]; - bank_array[i].xyz[1] = p->coord[0].xyz[1]; - bank_array[i].xyz[2] = p->coord[0].xyz[2]; + bank_array[i].r = p->r(); // Set that the bank particle is a neutron - bank_array[i].particle = static_cast(ParticleType::neutron); + bank_array[i].particle = Particle::Type::neutron; // Set the weight of the fission bank site bank_array[i].wgt = 1. / weight; @@ -168,23 +166,23 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, // Sample the azimuthal angle uniformly in [0, 2.pi) double phi = 2. * PI * prn(); - bank_array[i].uvw[0] = mu; - bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi); - bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi); + bank_array[i].u.x = mu; + bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi); + bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi); // Sample secondary energy distribution for the fission reaction and set // the energy in the fission bank int dg; int gout; - data::macro_xs[p->material].sample_fission_energy(p->g - 1, dg, gout); - bank_array[i].E = static_cast(gout + 1); + data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout); + bank_array[i].E = gout + 1; bank_array[i].delayed_group = dg + 1; // Set the delayed group on the particle as well - p->delayed_group = dg + 1; + p->delayed_group_ = dg + 1; // Increment the number of neutrons born delayed - if (p->delayed_group > 0) { + if (p->delayed_group_ > 0) { nu_d[dg]++; } } @@ -193,10 +191,10 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, *size_bank = std::min(*size_bank + nu, bank_array_size); // Store the total weight banked for analog fission tallies - p->n_bank = nu; - p->wgt_bank = nu / weight; + p->n_bank_ = nu; + p->wgt_bank_ = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p->n_delayed_bank[d] = nu_d[d]; + p->n_delayed_bank_[d] = nu_d[d]; } } @@ -205,26 +203,26 @@ absorption(Particle* p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * + p->wgt_absorb_ = p->wgt_ * simulation::material_xs.absorption / simulation::material_xs.total; // Adjust weight of particle by the probability of absorption - p->wgt -= p->absorb_wgt; - p->last_wgt = p->wgt; + p->wgt_ -= p->wgt_absorb_; + p->wgt_last_ = p->wgt_; // Score implicit absorpion estimate of keff #pragma omp atomic - global_tally_absorption += p->absorb_wgt * + global_tally_absorption += p->wgt_absorb_ * simulation::material_xs.nu_fission / simulation::material_xs.absorption; } else { if (simulation::material_xs.absorption > prn() * simulation::material_xs.total) { #pragma omp atomic - global_tally_absorption += p->wgt * simulation::material_xs.nu_fission / + global_tally_absorption += p->wgt_ * simulation::material_xs.nu_fission / simulation::material_xs.absorption; - p->alive = false; - p->event = EVENT_ABSORB; + p->alive_ = false; + p->event_ = EVENT_ABSORB; } } diff --git a/src/plot.cpp b/src/plot.cpp index ece9c76e7..ae2f00023 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -106,49 +106,48 @@ void create_ppm(Plot pl) data.resize({width, height}); int in_i, out_i; - double xyz[3]; + Position r; switch(pl.basis_) { case PlotBasis::xy : in_i = 0; out_i = 1; - xyz[0] = pl.origin_[0] - pl.width_[0] / 2.; - xyz[1] = pl.origin_[1] + pl.width_[1] / 2.; - xyz[2] = pl.origin_[2]; + r.x = pl.origin_[0] - pl.width_[0] / 2.; + r.y = pl.origin_[1] + pl.width_[1] / 2.; + r.z = pl.origin_[2]; break; case PlotBasis::xz : in_i = 0; out_i = 2; - xyz[0] = pl.origin_[0] - pl.width_[0] / 2.; - xyz[1] = pl.origin_[1]; - xyz[2] = pl.origin_[2] + pl.width_[1] / 2.; + r.x = pl.origin_[0] - pl.width_[0] / 2.; + r.y = pl.origin_[1]; + r.z = pl.origin_[2] + pl.width_[1] / 2.; break; case PlotBasis::yz : in_i = 1; out_i = 2; - xyz[0] = pl.origin_[0]; - xyz[1] = pl.origin_[1] - pl.width_[0] / 2.; - xyz[2] = pl.origin_[2] + pl.width_[1] / 2.; + r.x = pl.origin_[0]; + r.y = pl.origin_[1] - pl.width_[0] / 2.; + r.z = pl.origin_[2] + pl.width_[1] / 2.; break; } - double dir[3] = {0.5, 0.5, 0.5}; + Direction u {0.5, 0.5, 0.5}; #pragma omp parallel { Particle p; - p.initialize(); - std::copy(xyz, xyz+3, p.coord[0].xyz); - std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = model::root_universe; + p.r() = r; + p.u() = u; + p.coord_[0].universe = model::root_universe; #pragma omp for for (int y = 0; y < height; y++) { - p.coord[0].xyz[out_i] = xyz[out_i] - out_pixel * y; + p.r()[out_i] = r[out_i] - out_pixel * y; for (int x = 0; x < width; x++) { // local variables RGBColor rgb; int id; - p.coord[0].xyz[in_i] = xyz[in_i] + in_pixel * x; + p.r()[in_i] = r[in_i] + in_pixel * x; position_rgb(p, pl, rgb, id); data(x,y) = rgb; } @@ -309,11 +308,9 @@ void Plot::set_origin(pugi::xml_node plot_node) { // Copy plotting origin - std::vector pl_origin = get_node_array(plot_node, "origin"); + auto pl_origin = get_node_array(plot_node, "origin"); if (pl_origin.size() == 3) { - origin_[0] = pl_origin[0]; - origin_[1] = pl_origin[1]; - origin_[2] = pl_origin[2]; + origin_ = pl_origin; } else { std::stringstream err_msg; err_msg << "Origin must be length 3 in plot " @@ -329,8 +326,8 @@ Plot::set_width(pugi::xml_node plot_node) std::vector pl_width = get_node_array(plot_node, "width"); if (PlotType::slice == type_) { if (pl_width.size() == 2) { - width_[0] = pl_width[0]; - width_[1] = pl_width[1]; + width_.x = pl_width[0]; + width_.y = pl_width[1]; } else { std::stringstream err_msg; err_msg << " must be length 2 in slice plot " @@ -340,9 +337,7 @@ Plot::set_width(pugi::xml_node plot_node) } else if (PlotType::voxel == type_) { if (pl_width.size() == 3) { pl_width = get_node_array(plot_node, "width"); - width_[0] = pl_width[0]; - width_[1] = pl_width[1]; - width_[2] = pl_width[2]; + width_ = pl_width; } else { std::stringstream err_msg; err_msg << " must be length 3 in voxel plot " @@ -648,11 +643,11 @@ index_meshlines_mesh_(-1) void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) { - p.n_coord = 1; + p.n_coord_ = 1; bool found_cell = find_cell(&p, 0); - int j = p.n_coord - 1; + int j = p.n_coord_ - 1; if (settings::check_overlaps) {check_cell_overlap(&p);} @@ -666,23 +661,23 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) } else { if (PlotColorBy::mats == pl.color_by_) { // Assign color based on material - Cell* c = model::cells[p.coord[j].cell]; + const auto& c = model::cells[p.coord_[j].cell]; if (c->type_ == FILL_UNIVERSE) { // If we stopped on a middle universe level, treat as if not found rgb = pl.not_found_; id = -1; - } else if (p.material == MATERIAL_VOID) { + } else if (p.material_ == MATERIAL_VOID) { // By default, color void cells white rgb = WHITE; id = -1; } else { - rgb = pl.colors_[p.material]; - id = model::materials[p.material]->id_; + rgb = pl.colors_[p.material_]; + id = model::materials[p.material_]->id_; } } else if (PlotColorBy::cells == pl.color_by_) { // Assign color based on cell - rgb = pl.colors_[p.coord[j].cell]; - id = model::cells[p.coord[j].cell]->id_; + rgb = pl.colors_[p.coord_[j].cell]; + id = model::cells[p.coord_[j].cell]->id_; } } // endif found_cell } @@ -746,36 +741,27 @@ void draw_mesh_lines(Plot pl, ImageData& data) break; } - double xyz_ll_plot[3], xyz_ur_plot[3]; - xyz_ll_plot[0] = pl.origin_[0]; - xyz_ll_plot[1] = pl.origin_[1]; - xyz_ll_plot[2] = pl.origin_[2]; + Position ll_plot {pl.origin_}; + Position ur_plot {pl.origin_}; - xyz_ur_plot[0] = pl.origin_[0]; - xyz_ur_plot[1] = pl.origin_[1]; - xyz_ur_plot[2] = pl.origin_[2]; + ll_plot[outer] -= pl.width_[0] / 2.; + ll_plot[inner] -= pl.width_[1] / 2.; + ur_plot[outer] += pl.width_[0] / 2.; + ur_plot[inner] += pl.width_[1] / 2.; - xyz_ll_plot[outer] = pl.origin_[outer] - pl.width_[0] / 2.; - xyz_ll_plot[inner] = pl.origin_[inner] - pl.width_[1] / 2.; - xyz_ur_plot[outer] = pl.origin_[outer] + pl.width_[0] / 2.; - xyz_ur_plot[inner] = pl.origin_[inner] + pl.width_[1] / 2.; - - int width[3]; - width[0] = xyz_ur_plot[0] - xyz_ll_plot[0]; - width[1] = xyz_ur_plot[1] - xyz_ll_plot[1]; - width[2] = xyz_ur_plot[2] - xyz_ll_plot[2]; + Position width = ur_plot - ll_plot; auto& m = model::meshes[pl.index_meshlines_mesh_]; int ijk_ll[3], ijk_ur[3]; bool in_mesh; - m->get_indices(Position(xyz_ll_plot), &(ijk_ll[0]), &in_mesh); - m->get_indices(Position(xyz_ur_plot), &(ijk_ur[0]), &in_mesh); + m->get_indices(ll_plot, &(ijk_ll[0]), &in_mesh); + m->get_indices(ur_plot, &(ijk_ur[0]), &in_mesh); // Fortran/C++ index correction ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++; - double xyz_ll[3], xyz_ur[3]; + Position r_ll, r_ur; // sweep through all meshbins on this plane and draw borders for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) { for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) { @@ -783,20 +769,20 @@ void draw_mesh_lines(Plot pl, ImageData& data) if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) { int outrange[3], inrange[3]; // get xyz's of lower left and upper right of this mesh cell - xyz_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1); - xyz_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1); - xyz_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i; - xyz_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j; + r_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1); + r_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1); + r_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i; + r_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j; // map the xyz ranges to pixel ranges - double frac = (xyz_ll[outer] - xyz_ll_plot[outer]) / width[outer]; + double frac = (r_ll[outer] - ll_plot[outer]) / width[outer]; outrange[0] = int(frac * double(pl.pixels_[0])); - frac = (xyz_ur[outer] - xyz_ll_plot[outer]) / width[outer]; + frac = (r_ur[outer] - ll_plot[outer]) / width[outer]; outrange[1] = int(frac * double(pl.pixels_[0])); - frac = (xyz_ur[inner] - xyz_ll_plot[inner]) / width[inner]; + frac = (r_ur[inner] - ll_plot[inner]) / width[inner]; inrange[0] = int((1. - frac) * (double)pl.pixels_[1]); - frac = (xyz_ll[inner] - xyz_ll_plot[inner]) / width[inner]; + frac = (r_ll[inner] - ll_plot[inner]) / width[inner]; inrange[1] = int((1. - frac) * (double)pl.pixels_[1]); // draw lines @@ -846,18 +832,14 @@ void create_voxel(Plot pl) vox[2] = pl.width_[2]/(double)pl.pixels_[2]; // initial particle position - std::array ll; - ll[0] = pl.origin_[0] - pl.width_[0] / 2.; - ll[1] = pl.origin_[1] - pl.width_[1] / 2.; - ll[2] = pl.origin_[2] - pl.width_[2] / 2.; + Position ll = pl.origin_ - pl.width_ / 2.; // allocate and initialize particle - double dir[3] = {0.5, 0.5, 0.5}; + Direction u {0.5, 0.5, 0.5}; Particle p; - p.initialize(); - std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz); - std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = model::root_universe; + p.r() = ll; + p.u() = u; + p.coord_[0].universe = model::root_universe; // Open binary plot file for writing std::ofstream of; @@ -891,9 +873,9 @@ void create_voxel(Plot pl) voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace); // move to center of voxels - ll[0] = ll[0] + vox[0] / 2.; - ll[1] = ll[1] + vox[1] / 2.; - ll[2] = ll[2] + vox[2] / 2.; + ll.x += vox[0] / 2.; + ll.y += vox[1] / 2.; + ll.z += vox[2] / 2.; int data[pl.pixels_[1]][pl.pixels_[0]]; @@ -910,16 +892,16 @@ void create_voxel(Plot pl) // write to plot data data[y][x] = id; // advance particle in x direction - p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0]; + p.r().x += vox[0]; } // advance particle in y direction - p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1]; - p.coord[0].xyz[0] = ll[0]; + p.r().y += vox[1]; + p.r().x = ll[0]; } // advance particle in z direction - p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2]; - p.coord[0].xyz[1] = ll[1]; - p.coord[0].xyz[0] = ll[0]; + p.r().z += vox[2]; + p.r().y = ll[1]; + p.r().x = ll[0]; // Write to HDF5 dataset voxel_write_slice(z, dspace, dset, memspace, &(data[0])); } diff --git a/src/position.cpp b/src/position.cpp index ecf5be631..8ab97c96d 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -78,4 +78,10 @@ Position::operator/=(double v) return *this; } +Position +Position::operator-() const +{ + return {-x, -y, -z}; +} + } // namespace openmc diff --git a/src/reaction.cpp b/src/reaction.cpp index 0fcac2e1d..1f997f078 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -72,7 +72,7 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) // mark fission reactions so that we avoid the angle sampling. if (is_fission(mt_)) { for (auto& p : products_) { - if (p.particle_ == ParticleType::neutron) { + if (p.particle_ == Particle::Type::neutron) { for (auto& d : p.distribution_) { auto d_ = dynamic_cast(d.get()); if (d_) d_->fission() = true; diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 7e3175343..4f3994ff2 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -22,9 +22,9 @@ ReactionProduct::ReactionProduct(hid_t group) std::string temp; read_attribute(group, "particle", temp); if (temp == "neutron") { - particle_ = ParticleType::neutron; + particle_ = Particle::Type::neutron; } else if (temp == "photon") { - particle_ = ParticleType::photon; + particle_ = Particle::Type::photon; } // Read emission mode and decay rate @@ -69,13 +69,13 @@ ReactionProduct::ReactionProduct(hid_t group) // Determine distribution type and read data read_attribute(dgroup, "type", temp); if (temp == "uncorrelated") { - distribution_.emplace_back(new UncorrelatedAngleEnergy{dgroup}); + distribution_.push_back(std::make_unique(dgroup)); } else if (temp == "correlated") { - distribution_.emplace_back(new CorrelatedAngleEnergy{dgroup}); + distribution_.push_back(std::make_unique(dgroup)); } else if (temp == "nbody") { - distribution_.emplace_back(new NBodyPhaseSpace{dgroup}); + distribution_.push_back(std::make_unique(dgroup)); } else if (temp == "kalbach-mann") { - distribution_.emplace_back(new KalbachMann{dgroup}); + distribution_.push_back(std::make_unique(dgroup)); } close_group(dgroup); diff --git a/src/settings.cpp b/src/settings.cpp index 92e85ab17..dfa5ebb7e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -509,7 +509,7 @@ void read_settings_xml() // Read entropy mesh from auto node_entropy = root.child("entropy"); - model::meshes.emplace_back(new RegularMesh{node_entropy}); + model::meshes.push_back(std::make_unique(node_entropy)); // Set entropy mesh index index_entropy_mesh = model::meshes.size() - 1; @@ -555,7 +555,7 @@ void read_settings_xml() // Read entropy mesh from auto node_ufs = root.child("uniform_fs"); - model::meshes.emplace_back(new RegularMesh{node_ufs}); + model::meshes.push_back(std::make_unique(node_ufs)); // Set entropy mesh index index_ufs_mesh = model::meshes.size() - 1; diff --git a/src/simulation.cpp b/src/simulation.cpp index 50ce1b5c0..7d30ccea9 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -474,29 +474,29 @@ void initialize_history(Particle* p, int64_t index_source) p->from_source(&simulation::source_bank[index_source - 1]); // set identifier for particle - p->id = simulation::work_index[mpi::rank] + index_source; + p->id_ = simulation::work_index[mpi::rank] + index_source; // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) - * settings::n_particles + p->id; + * settings::n_particles + p->id_; set_particle_seed(particle_seed); // set particle trace simulation::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) simulation::trace = true; // Set particle track. - p->write_track = false; + p->write_track_ = false; if (settings::write_all_tracks) { - p->write_track = true; + p->write_track_ = true; } else if (settings::track_identifiers.size() > 0) { for (const auto& t : settings::track_identifiers) { if (simulation::current_batch == t[0] && simulation::current_gen == t[1] && - p->id == t[2]) { - p->write_track = true; + p->id_ == t[2]) { + p->write_track_ = true; break; } } diff --git a/src/source.cpp b/src/source.cpp index 1ebdc235b..36d66394c 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -47,9 +47,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) if (check_for_node(node, "particle")) { auto temp_str = get_node_value(node, "particle", true, true); if (temp_str == "neutron") { - particle_ = ParticleType::neutron; + particle_ = Particle::Type::neutron; } else if (temp_str == "photon") { - particle_ = ParticleType::photon; + particle_ = Particle::Type::photon; settings::photon_transport = true; } else { fatal_error(std::string("Unknown source particle type: ") + temp_str); @@ -140,9 +140,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) } -Bank SourceDistribution::sample() const +Particle::Bank SourceDistribution::sample() const { - Bank site; + Particle::Bank site; // Set weight to one by default site.wgt = 1.0; @@ -153,17 +153,15 @@ Bank SourceDistribution::sample() const static int n_accept = 0; while (!found) { // Set particle type - site.particle = static_cast(particle_); + site.particle = particle_; // Sample spatial distribution - Position r = space_->sample(); - site.xyz[0] = r.x; - site.xyz[1] = r.y; - site.xyz[2] = r.z; + site.r = space_->sample(); + double xyz[] {site.r.x, site.r.y, site.r.z}; // Now search to see if location exists in geometry int32_t cell_index, instance; - int err = openmc_find_cell(site.xyz, &cell_index, &instance); + int err = openmc_find_cell(xyz, &cell_index, &instance); found = (err != OPENMC_E_GEOMETRY); // Check if spatial site is in fissionable material @@ -172,7 +170,7 @@ Bank SourceDistribution::sample() const if (space_box) { if (space_box->only_fissionable()) { // Determine material - auto c = model::cells[cell_index]; + const auto& c = model::cells[cell_index]; auto mat_index = c->material_.size() == 1 ? c->material_[0] : c->material_[instance]; @@ -200,10 +198,7 @@ Bank SourceDistribution::sample() const ++n_accept; // Sample angle - Direction u = angle_->sample(); - site.uvw[0] = u.x; - site.uvw[1] = u.y; - site.uvw[2] = u.z; + site.u = angle_->sample(); // Check for monoenergetic source above maximum particle energy auto p = static_cast(particle_); @@ -290,7 +285,7 @@ void initialize_source() } } -Bank sample_external_source() +Particle::Bank sample_external_source() { // Set the random number generator to the source stream. prn_set_stream(STREAM_SOURCE); @@ -312,7 +307,7 @@ Bank sample_external_source() } // Sample source site from i-th source distribution - Bank site {model::external_sources[i].sample()}; + Particle::Bank site {model::external_sources[i].sample()}; // If running in MG, convert site % E to group if (!settings::run_CE) { diff --git a/src/state_point.cpp b/src/state_point.cpp index cca622950..ab1473403 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -467,19 +467,21 @@ void load_state_point() hid_t h5banktype() { - // Create type for array of 3 reals - hsize_t dims[] {3}; - hid_t triplet = H5Tarray_create(H5T_NATIVE_DOUBLE, 1, dims); + // Create compound type for position + hid_t postype = H5Tcreate(H5T_COMPOUND, sizeof(struct Position)); + H5Tinsert(postype, "x", HOFFSET(Position, x), H5T_NATIVE_DOUBLE); + H5Tinsert(postype, "y", HOFFSET(Position, y), H5T_NATIVE_DOUBLE); + H5Tinsert(postype, "z", HOFFSET(Position, z), H5T_NATIVE_DOUBLE); // Create bank datatype - hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Bank)); - H5Tinsert(banktype, "wgt", HOFFSET(Bank, wgt), H5T_NATIVE_DOUBLE); - H5Tinsert(banktype, "xyz", HOFFSET(Bank, xyz), triplet); - H5Tinsert(banktype, "uvw", HOFFSET(Bank, uvw), triplet); - H5Tinsert(banktype, "E", HOFFSET(Bank, E), H5T_NATIVE_DOUBLE); - H5Tinsert(banktype, "delayed_group", HOFFSET(Bank, delayed_group), H5T_NATIVE_INT); + hid_t banktype = H5Tcreate(H5T_COMPOUND, sizeof(struct Particle::Bank)); + H5Tinsert(banktype, "r", HOFFSET(Particle::Bank, r), postype); + H5Tinsert(banktype, "u", HOFFSET(Particle::Bank, u), postype); + H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE); + H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT); - H5Tclose(triplet); + H5Tclose(postype); return banktype; } @@ -565,7 +567,7 @@ write_source_bank(hid_t group_id) // Save source bank sites since the souce_bank array is overwritten below #ifdef OPENMC_MPI - std::vector temp_source {simulation::source_bank.begin(), + std::vector temp_source {simulation::source_bank.begin(), simulation::source_bank.begin() + simulation::work}; #endif diff --git a/src/summary.cpp b/src/summary.cpp index 1eeab5d04..500a5c15c 100644 --- a/src/summary.cpp +++ b/src/summary.cpp @@ -102,19 +102,19 @@ void write_geometry(hid_t file) write_attribute(geom_group, "n_lattices", model::lattices.size()); auto cells_group = create_group(geom_group, "cells"); - for (Cell* c : model::cells) c->to_hdf5(cells_group); + for (const auto& c : model::cells) c->to_hdf5(cells_group); close_group(cells_group); auto surfaces_group = create_group(geom_group, "surfaces"); - for (Surface* surf : model::surfaces) surf->to_hdf5(surfaces_group); + for (const auto& surf : model::surfaces) surf->to_hdf5(surfaces_group); close_group(surfaces_group); auto universes_group = create_group(geom_group, "universes"); - for (Universe* u : model::universes) u->to_hdf5(universes_group); + for (const auto& u : model::universes) u->to_hdf5(universes_group); close_group(universes_group); auto lattices_group = create_group(geom_group, "lattices"); - for (Lattice* lat : model::lattices) lat->to_hdf5(lattices_group); + for (const auto& lat : model::lattices) lat->to_hdf5(lattices_group); close_group(lattices_group); close_group(geom_group); diff --git a/src/surface.cpp b/src/surface.cpp index d56adbfd9..bf0315698 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -27,10 +27,8 @@ extern "C" const int BC_PERIODIC {3}; //============================================================================== namespace model { - -std::vector surfaces; -std::map surface_map; - + std::vector> surfaces; + std::unordered_map surface_map; } // namespace model //============================================================================== @@ -1080,40 +1078,40 @@ void read_surfaces(pugi::xml_node node) std::string surf_type = get_node_value(surf_node, "type", true, true); if (surf_type == "x-plane") { - model::surfaces.push_back(new SurfaceXPlane(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "y-plane") { - model::surfaces.push_back(new SurfaceYPlane(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "z-plane") { - model::surfaces.push_back(new SurfaceZPlane(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "plane") { - model::surfaces.push_back(new SurfacePlane(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "x-cylinder") { - model::surfaces.push_back(new SurfaceXCylinder(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "y-cylinder") { - model::surfaces.push_back(new SurfaceYCylinder(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "z-cylinder") { - model::surfaces.push_back(new SurfaceZCylinder(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "sphere") { - model::surfaces.push_back(new SurfaceSphere(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "x-cone") { - model::surfaces.push_back(new SurfaceXCone(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "y-cone") { - model::surfaces.push_back(new SurfaceYCone(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "z-cone") { - model::surfaces.push_back(new SurfaceZCone(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else if (surf_type == "quadric") { - model::surfaces.push_back(new SurfaceQuadric(surf_node)); + model::surfaces.push_back(std::make_unique(surf_node)); } else { std::stringstream err_msg; @@ -1143,8 +1141,8 @@ void read_surfaces(pugi::xml_node node) for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) { if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = model::surfaces[i_surf]; - PeriodicSurface* surf = dynamic_cast(surf_base); + Surface* surf_base = model::surfaces[i_surf].get(); + auto surf = dynamic_cast(surf_base); // Make sure this surface inherits from PeriodicSurface. if (!surf) { @@ -1188,8 +1186,8 @@ void read_surfaces(pugi::xml_node node) for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) { if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) { // Downcast to the PeriodicSurface type. - Surface* surf_base = model::surfaces[i_surf]; - PeriodicSurface* surf = dynamic_cast(surf_base); + Surface* surf_base = model::surfaces[i_surf].get(); + auto surf = dynamic_cast(surf_base); // Also try downcasting to the SurfacePlane type (which must be handled // differently). @@ -1260,7 +1258,6 @@ void read_surfaces(pugi::xml_node node) void free_memory_surfaces() { - for (Surface* surf : model::surfaces) {delete surf;} model::surfaces.clear(); model::surface_map.clear(); } diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 2ad251b73..315cc081e 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -126,11 +126,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, if (score_bin == SCORE_FLUX) { score *= flux_deriv; return; - } else if (p->material == MATERIAL_VOID) { + } else if (p->material_ == MATERIAL_VOID) { score *= flux_deriv; return; } - const Material& material {*model::materials[p->material]}; + const Material& material {*model::materials[p->material_]}; if (material.id_ != deriv.diff_material) { score *= flux_deriv; return; @@ -190,7 +190,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (tally.estimator_) { case ESTIMATOR_ANALOG: - if (p->event_nuclide != deriv.diff_nuclide) { + if (p->event_nuclide_ != deriv.diff_nuclide) { score *= flux_deriv; return; } @@ -320,10 +320,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, // Find the index of the event nuclide. int i; for (i = 0; i < material.nuclide_.size(); ++i) - if (material.nuclide_[i] == p->event_nuclide) break; + if (material.nuclide_[i] == p->event_nuclide_) break; - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (!multipole_in_range(&nuc, p->last_E)) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (!multipole_in_range(&nuc, p->E_last_)) { score *= flux_deriv; break; } @@ -331,10 +331,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (simulation::micro_xs[p->event_nuclide].total) { + if (simulation::micro_xs[p->event_nuclide_].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) / simulation::material_xs.total; } else { @@ -343,11 +343,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (simulation::micro_xs[p->event_nuclide].total - - simulation::micro_xs[p->event_nuclide].absorption) { + if (simulation::micro_xs[p->event_nuclide_].total + - simulation::micro_xs[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_s * material.atom_density_(i) / (simulation::material_xs.total - simulation::material_xs.absorption); @@ -357,10 +357,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (simulation::micro_xs[p->event_nuclide].absorption) { + if (simulation::micro_xs[p->event_nuclide_].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_a * material.atom_density_(i) / simulation::material_xs.absorption; } else { @@ -369,10 +369,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (simulation::micro_xs[p->event_nuclide].fission) { + if (simulation::micro_xs[p->event_nuclide_].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_f * material.atom_density_(i) / simulation::material_xs.fission; } else { @@ -381,12 +381,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_NU_FISSION: - if (simulation::micro_xs[p->event_nuclide].fission) { - double nu = simulation::micro_xs[p->event_nuclide].nu_fission - / simulation::micro_xs[p->event_nuclide].fission; + if (simulation::micro_xs[p->event_nuclide_].fission) { + double nu = simulation::micro_xs[p->event_nuclide_].nu_fission + / simulation::micro_xs[p->event_nuclide_].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + nu * dsig_f * material.atom_density_(i) / simulation::material_xs.nu_fission; } else { @@ -404,7 +404,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case ESTIMATOR_COLLISION: if (i_nuclide != -1) { const auto& nuc {data::nuclides[i_nuclide]}; - if (!multipole_in_range(nuc.get(), p->last_E)) { + if (!multipole_in_range(nuc.get(), p->E_last_)) { score *= flux_deriv; return; } @@ -418,11 +418,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(&nuc, p->last_E) + if (multipole_in_range(&nuc, p->E_last_) && simulation::micro_xs[i_nuc].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); } } @@ -431,7 +431,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total; } else { @@ -446,12 +446,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(&nuc, p->last_E) + if (multipole_in_range(&nuc, p->E_last_) && (simulation::micro_xs[i_nuc].total - simulation::micro_xs[i_nuc].absorption)) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_s * material.atom_density_(i); } } @@ -462,7 +462,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total - simulation::micro_xs[i_nuclide].absorption); } else { @@ -476,11 +476,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(&nuc, p->last_E) + if (multipole_in_range(&nuc, p->E_last_) && simulation::micro_xs[i_nuc].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_a * material.atom_density_(i); } } @@ -489,7 +489,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_a / simulation::micro_xs[i_nuclide].absorption; } else { @@ -503,11 +503,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(&nuc, p->last_E) + if (multipole_in_range(&nuc, p->E_last_) && simulation::micro_xs[i_nuc].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += dsig_f * material.atom_density_(i); } } @@ -516,7 +516,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_f / simulation::micro_xs[i_nuclide].fission; } else { @@ -530,13 +530,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(&nuc, p->last_E) + if (multipole_in_range(&nuc, p->E_last_) && simulation::micro_xs[i_nuc].fission) { double nu = simulation::micro_xs[i_nuc].nu_fission / simulation::micro_xs[i_nuc].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); cum_dsig += nu * dsig_f * material.atom_density_(i); } } @@ -545,7 +545,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); score *= flux_deriv + dsig_f / simulation::micro_xs[i_nuclide].fission; } else { @@ -570,8 +570,8 @@ void score_track_derivative(const 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]}; + if (p->material_ == MATERIAL_VOID) return; + const Material& material {*model::materials[p->material_]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; @@ -597,13 +597,13 @@ score_track_derivative(const Particle* p, double distance) case DIFF_TEMPERATURE: for (auto i = 0; i < material.nuclide_.size(); ++i) { const auto& nuc {*data::nuclides[material.nuclide_[i]]}; - if (multipole_in_range(&nuc, p->last_E)) { + if (multipole_in_range(&nuc, p->E_last_)) { // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist // (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_); deriv.flux_deriv -= distance * (dsig_s + dsig_a) * material.atom_density_(i); } @@ -616,8 +616,8 @@ score_track_derivative(const Particle* p, double distance) void score_collision_derivative(const 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]}; + if (p->material_ == MATERIAL_VOID) return; + const Material& material {*model::materials[p->material_]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; @@ -632,7 +632,7 @@ void score_collision_derivative(const Particle* p) break; case DIFF_NUCLIDE_DENSITY: - if (p->event_nuclide != deriv.diff_nuclide) continue; + if (p->event_nuclide_ != deriv.diff_nuclide) continue; // Find the index in this material for the diff_nuclide. int i; for (i = 0; i < material.nuclide_.size(); ++i) @@ -656,19 +656,19 @@ void score_collision_derivative(const Particle* p) // Loop over the material's nuclides until we find the event nuclide. for (auto i_nuc : material.nuclide_) { const auto& nuc {*data::nuclides[i_nuc]}; - if (i_nuc == p->event_nuclide && multipole_in_range(&nuc, p->last_E)) { + if (i_nuc == p->event_nuclide_ && multipole_in_range(&nuc, p->E_last_)) { // phi is proportional to Sigma_s // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s const auto& micro_xs {simulation::micro_xs[i_nuc]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_); deriv.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, uvw'->uvw) = Sigma_s(E') * P(E'->E, uvw'->uvw). - // We are assuming that d_P(E'->E, uvw'->uvw) / d_T = 0 and only + // Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u). + // We are assuming that d_P(E'->E, u'->u) / d_T = 0 and only // computing d_S(E') / d_T. Using this approximation in the vicinity // of low-energy resonances causes errors (~2-5% for PWR pincell // eigenvalue derivatives). diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp index f4e265dd8..1a33a6bfe 100644 --- a/src/tallies/filter_azimuthal.cpp +++ b/src/tallies/filter_azimuthal.cpp @@ -42,9 +42,9 @@ AzimuthalFilter::get_all_bins(const Particle* p, int estimator, { double phi; if (estimator == ESTIMATOR_TRACKLENGTH) { - phi = std::atan2(p->coord[0].uvw[1], p->coord[0].uvw[0]); + phi = std::atan2(p->u().y, p->u().x); } else { - phi = std::atan2(p->last_uvw[1], p->last_uvw[0]); + phi = std::atan2(p->u_last_.y, p->u_last_.x); } if (phi >= bins_.front() && phi <= bins_.back()) { diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 1b2b05416..8286a6976 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -41,8 +41,8 @@ void CellFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - for (int i = 0; i < p->n_coord; i++) { - auto search = map_.find(p->coord[i].cell); + for (int i = 0; i < p->n_coord_; i++) { + auto search = map_.find(p->coord_[i].cell); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp index 005d3844b..150962297 100644 --- a/src/tallies/filter_cellborn.cpp +++ b/src/tallies/filter_cellborn.cpp @@ -8,7 +8,7 @@ void CellbornFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - auto search = map_.find(p->cell_born); + auto search = map_.find(p->cell_born_); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp index 9a43f7c43..b8e5d0b07 100644 --- a/src/tallies/filter_cellfrom.cpp +++ b/src/tallies/filter_cellfrom.cpp @@ -8,8 +8,8 @@ void CellFromFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - for (int i = 0; i < p->last_n_coord; i++) { - auto search = map_.find(p->last_cell[i]); + for (int i = 0; i < p->n_coord_last_; i++) { + auto search = map_.find(p->cell_last_[i]); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index ab7ea67cf..94f23621b 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -40,20 +40,20 @@ DistribcellFilter::get_all_bins(const Particle* p, int estimator, { int offset = 0; auto distribcell_index = model::cells[cell_]->distribcell_index_; - for (int i = 0; i < p->n_coord; i++) { - auto& c {*model::cells[p->coord[i].cell]}; + for (int i = 0; i < p->n_coord_; i++) { + auto& c {*model::cells[p->coord_[i].cell]}; if (c.type_ == FILL_UNIVERSE) { offset += c.offset_[distribcell_index]; } else if (c.type_ == FILL_LATTICE) { - auto& lat {*model::lattices[p->coord[i+1].lattice-1]}; - int i_xyz[3] {p->coord[i+1].lattice_x, - p->coord[i+1].lattice_y, - p->coord[i+1].lattice_z}; + auto& lat {*model::lattices[p->coord_[i+1].lattice]}; + int i_xyz[3] {p->coord_[i+1].lattice_x, + p->coord_[i+1].lattice_y, + p->coord_[i+1].lattice_z}; if (lat.are_valid_indices(i_xyz)) { offset += lat.offset(distribcell_index, i_xyz); } } - if (cell_ == p->coord[i].cell) { + if (cell_ == p->coord_[i].cell) { match.bins_.push_back(offset); match.weights_.push_back(1.0); return; diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index 4cbde1cc4..e11ab3080 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -41,17 +41,17 @@ void EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->g != F90_NONE && matches_transport_groups_) { + if (p->g_ != F90_NONE && matches_transport_groups_) { if (estimator == ESTIMATOR_TRACKLENGTH) { - match.bins_.push_back(data::num_energy_groups - p->g); + match.bins_.push_back(data::num_energy_groups - p->g_); } else { - match.bins_.push_back(data::num_energy_groups - p->last_g); + match.bins_.push_back(data::num_energy_groups - p->g_last_); } match.weights_.push_back(1.0); } else { // Get the pre-collision energy of the particle. - auto E = p->last_E; + auto E = p->E_last_; // Bin the energy. if (E >= bins_.front() && E <= bins_.back()) { @@ -85,13 +85,13 @@ void EnergyoutFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->g != F90_NONE && matches_transport_groups_) { - match.bins_.push_back(data::num_energy_groups - p->g); + if (p->g_ != F90_NONE && matches_transport_groups_) { + match.bins_.push_back(data::num_energy_groups - p->g_); match.weights_.push_back(1.0); } else { - if (p->E >= bins_.front() && p->E <= bins_.back()) { - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E); + if (p->E_ >= bins_.front() && p->E_ <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E_); match.bins_.push_back(bin); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 3c616aabf..9b28e0f76 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -33,12 +33,12 @@ void EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) { + if (p->E_last_ >= energy_.front() && p->E_last_ <= energy_.back()) { // Search for the incoming energy bin. - auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E); + auto i = lower_bound_index(energy_.begin(), energy_.end(), p->E_last_); // Compute the interpolation factor between the nearest bins. - double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]); + double f = (p->E_last_ - energy_[i]) / (energy_[i+1] - energy_[i]); // Interpolate on the lin-lin grid. match.bins_.push_back(0); diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp index 1e799932a..1d1d440b2 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -19,7 +19,7 @@ LegendreFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { double wgt[n_bins_]; - calc_pn_c(order_, p->mu, wgt); + calc_pn_c(order_, p->mu_, wgt); for (int i = 0; i < n_bins_; i++) { match.bins_.push_back(i); match.weights_.push_back(wgt[i]); diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index a4da2deb1..fd9e637bd 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -42,7 +42,7 @@ void MaterialFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - auto search = map_.find(p->material); + auto search = map_.find(p->material_); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 177dc407e..f69a96652 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -35,7 +35,7 @@ MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { if (estimator != ESTIMATOR_TRACKLENGTH) { - auto bin = model::meshes[mesh_]->get_bin(p->coord[0].xyz); + auto bin = model::meshes[mesh_]->get_bin(p->r()); if (bin >= 0) { match.bins_.push_back(bin); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp index 78a648e23..aed8371b8 100644 --- a/src/tallies/filter_mu.cpp +++ b/src/tallies/filter_mu.cpp @@ -38,8 +38,8 @@ void MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->mu >= bins_.front() && p->mu <= bins_.back()) { - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu); + if (p->mu_ >= bins_.front() && p->mu_ <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu_); match.bins_.push_back(bin); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index 967febc06..a2902785a 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -7,8 +7,10 @@ namespace openmc { void ParticleFilter::from_xml(pugi::xml_node node) { - particles_ = get_node_array(node, "bins"); - for (auto& p : particles_) --p; + auto particles = get_node_array(node, "bins"); + for (auto& p : particles) { + particles_.push_back(static_cast(p - 1)); + } n_bins_ = particles_.size(); } @@ -17,7 +19,7 @@ ParticleFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { for (auto i = 0; i < particles_.size(); i++) { - if (particles_[i] == p->type) { + if (particles_[i] == p->type_) { match.bins_.push_back(i); match.weights_.push_back(1.0); } @@ -28,13 +30,26 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); - write_dataset(filter_group, "bins", particles_); + std::vector particles; + for (auto p : particles_) { + particles.push_back(static_cast(p) + 1); + } + write_dataset(filter_group, "bins", particles); } std::string ParticleFilter::text_label(int bin) const { - return "Particle " + std::to_string(particles_[bin]); + switch (particles_[bin]) { + case Particle::Type::neutron: + return "Particle: neutron"; + case Particle::Type::photon: + return "Particle: photon"; + case Particle::Type::electron: + return "Particle: electron"; + case Particle::Type::positron: + return "Particle: positron"; + } } } // namespace openmc diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp index 5216511a9..9730bbb58 100644 --- a/src/tallies/filter_polar.cpp +++ b/src/tallies/filter_polar.cpp @@ -41,9 +41,9 @@ const { double theta; if (estimator == ESTIMATOR_TRACKLENGTH) { - theta = std::acos(p->coord[0].uvw[2]); + theta = std::acos(p->u().z); } else { - theta = std::acos(p->last_uvw[2]); + theta = std::acos(p->u_last_.z); } if (theta >= bins_.front() && theta <= bins_.back()) { diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index 253bf7ff3..9aae06bee 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -37,14 +37,14 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator, // Determine cosine term for scatter expansion if necessary double wgt[order_ + 1]; if (cosine_ == SphericalHarmonicsCosine::scatter) { - calc_pn_c(order_, p->mu, wgt); + calc_pn_c(order_, p->mu_, wgt); } else { for (int i = 0; i < order_ + 1; i++) wgt[i] = 1; } // Find the Rn,m values double rn[n_bins_]; - calc_rn_c(order_, p->last_uvw, rn); + calc_rn(order_, p->u_last_, rn); int j = 0; for (int n = 0; n < order_ + 1; n++) { diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index dae62d1d2..616dbbbd4 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -38,11 +38,11 @@ SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator, // Get the coordinate along the axis of interest. double x; if (axis_ == LegendreAxis::x) { - x = p->coord[0].xyz[0]; + x = p->r().x; } else if (axis_ == LegendreAxis::y) { - x = p->coord[0].xyz[1]; + x = p->r().y; } else { - x = p->coord[0].xyz[2]; + x = p->r().z; } if (x >= min_ && x <= max_) { diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index 774571d75..faa8aafbf 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -41,10 +41,10 @@ void SurfaceFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - auto search = map_.find(std::abs(p->surface)-1); + auto search = map_.find(std::abs(p->surface_)-1); if (search != map_.end()) { match.bins_.push_back(search->second); - if (p->surface < 0) { + if (p->surface_ < 0) { match.weights_.push_back(-1.0); } else { match.weights_.push_back(1.0); diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp index 5a6d5e1eb..50c058c05 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -41,8 +41,8 @@ void UniverseFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - for (int i = 0; i < p->n_coord; i++) { - auto search = map_.find(p->coord[i].universe); + for (int i = 0; i < p->n_coord_; i++) { + auto search = map_.find(p->coord_[i].universe); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index a005f7ddc..c4e620c76 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -29,8 +29,8 @@ ZernikeFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { // Determine the normalized (r,theta) coordinates. - double x = p->coord[0].xyz[0] - x_; - double y = p->coord[0].xyz[1] - y_; + double x = p->r().x - x_; + double y = p->r().y - y_; double r = std::sqrt(x*x + y*y) / r_; double theta = std::atan2(y, x); @@ -86,8 +86,8 @@ ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { // Determine the normalized radius coordinate. - double x = p->coord[0].xyz[0] - x_; - double y = p->coord[0].xyz[1] - y_; + double x = p->r().x - x_; + double y = p->r().y - y_; double r = std::sqrt(x*x + y*y) / r_; if (r <= 1.0) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 4d55350a0..5bae26b75 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -765,9 +765,9 @@ void read_tallies_xml() } else { const auto& f = model::tally_filters[particle_filter_index].get(); auto pf = dynamic_cast(f); - for (int p : pf->particles_) { - if (p == static_cast(ParticleType::electron) || - p == static_cast(ParticleType::positron)) { + for (auto p : pf->particles_) { + if (p == Particle::Type::electron || + p == Particle::Type::positron) { t->estimator_ = ESTIMATOR_ANALOG; } } @@ -776,11 +776,11 @@ void read_tallies_xml() if (particle_filter_index >= 0) { const auto& f = model::tally_filters[particle_filter_index].get(); auto pf = dynamic_cast(f); - for (int p : pf->particles_) { - if (p != static_cast(ParticleType::neutron)) { + for (auto p : pf->particles_) { + if (p != Particle::Type::neutron) { warning("Particle filter other than NEUTRON used with photon " "transport turned off. All tallies for particle type " + - std::to_string(p) + " will have no scores"); + std::to_string(static_cast(p)) + " will have no scores"); } } } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5288ea597..ed4d09254 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1,5 +1,6 @@ #include "openmc/tallies/tally_scoring.h" +#include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -190,8 +191,8 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) // rate. Otherwise, the sum of all nu-fission rates would be ~1.0. // loop over number of particles banked - for (auto i = 0; i < p->n_bank; ++i) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; // get the delayed group @@ -322,7 +323,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto& tally {*model::tallies[i_tally]}; // Get the pre-collision energy of the particle. - auto E = p->last_E; + auto E = p->E_last_; for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; @@ -341,13 +342,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->wgt_last_ + p->wgt_absorb_; } else { - score = p->last_wgt; + score = p->wgt_last_; } - if (p->type == static_cast(ParticleType::neutron) || - p->type == static_cast(ParticleType::photon)) { + if (p->type_ == Particle::Type::neutron || + p->type_ == Particle::Type::photon) { score *= flux / simulation::material_xs.total; } else { score = 0.; @@ -365,9 +366,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = (p->last_wgt + p->absorb_wgt) * flux; + score = (p->wgt_last_ + p->wgt_absorb_) * flux; } else { - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } } else { @@ -388,9 +389,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->wgt_last_ + p->wgt_absorb_; } else { - score = p->last_wgt; + score = p->wgt_last_; } score *= flux / simulation::material_xs.total; } else { @@ -404,10 +405,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_SCATTER) continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { score = (simulation::micro_xs[i_nuclide].total @@ -423,21 +424,21 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_NU_SCATTER: // Only analog estimators are available. // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_SCATTER) continue; // For scattering production, we need to use the pre-collision weight // times the yield as the estimate for the number of neutrons exiting a // reaction with neutrons in the exit channel - if (p->event_MT == ELASTIC || p->event_MT == N_LEVEL || - (p->event_MT >= N_N1 && p->event_MT <= N_NC)) { + if (p->event_mt_ == ELASTIC || p->event_mt_ == N_LEVEL || + (p->event_mt_ >= N_N1 && p->event_mt_ <= N_NC)) { // Don't waste time on very common reactions we know have // multiplicities of one. - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } else { // Get yield and apply to score auto m = - data::nuclides[p->event_nuclide]->reaction_index_[p->event_MT]; - const auto& rxn {*data::nuclides[p->event_nuclide]->reactions_[m]}; - score = p->last_wgt * flux * (*rxn.products_[0].yield_)(E); + data::nuclides[p->event_nuclide_]->reaction_index_[p->event_mt_]; + const auto& rxn {*data::nuclides[p->event_nuclide_]->reactions_[m]}; + score = p->wgt_last_ * flux * (*rxn.products_[0].yield_)(E); } break; @@ -447,13 +448,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } } else { if (i_nuclide >= 0) { @@ -473,22 +474,22 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->wgt_absorb_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->last_wgt - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_last_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { if (i_nuclide >= 0) { @@ -503,7 +504,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -515,22 +516,22 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].nu_fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->wgt_absorb_ + * simulation::micro_xs[p->event_nuclide_].nu_fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p->wgt_bank * flux; + score = simulation::keff * p->wgt_bank_ * flux; } } else { if (i_nuclide >= 0) { @@ -546,7 +547,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_PROMPT_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -558,27 +559,27 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - * data::nuclides[p->event_nuclide] + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->wgt_absorb_ + * simulation::micro_xs[p->event_nuclide_].fission + * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::prompt) - / simulation::micro_xs[p->event_nuclide].absorption * flux; + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p->n_delayed_bank, - p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank); - score = simulation::keff * p->wgt_bank * prompt_frac * flux; + auto n_delayed = std::accumulate(p->n_delayed_bank_, + p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank_); + score = simulation::keff * p->wgt_bank_ * prompt_frac * flux; } } else { if (i_nuclide >= 0) { @@ -589,8 +590,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = 0.; // Add up contributions from each nuclide in the material. - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -608,7 +609,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -620,8 +621,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0 - && data::nuclides[p->event_nuclide]->fissionable_) { + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 + && data::nuclides[p->event_nuclide_]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -630,11 +631,11 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - auto yield = data::nuclides[p->event_nuclide] + auto yield = data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_absorb_ * yield + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -643,16 +644,16 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - * data::nuclides[p->event_nuclide] + score = p->wgt_absorb_ + * simulation::micro_xs[p->event_nuclide_].fission + * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed) - / simulation::micro_xs[p->event_nuclide].absorption *flux; + / simulation::micro_xs[p->event_nuclide_].absorption *flux; } } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -668,16 +669,16 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - score = simulation::keff * p->wgt_bank / p->n_bank - * p->n_delayed_bank[d-1] * flux; + 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); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p->n_delayed_bank, - p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p->wgt_bank / p->n_bank * n_delayed + auto n_delayed = std::accumulate(p->n_delayed_bank_, + p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + score = simulation::keff * p->wgt_bank_ / p->n_bank_ * n_delayed * flux; } } @@ -712,8 +713,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -732,8 +733,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, continue; } else { score = 0.; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -756,8 +757,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -771,9 +772,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption + score = p->wgt_absorb_ * yield + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -794,15 +795,15 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += rate * p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission * yield - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score += rate * p->wgt_absorb_ + * simulation::micro_xs[p->event_nuclide_].fission * yield + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -811,12 +812,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p->n_bank; ++i) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { - const auto& nuc {*data::nuclides[p->event_nuclide]}; + const auto& nuc {*data::nuclides[p->event_nuclide_]}; const auto& rxn {*nuc.fission_rx_[0]}; auto rate = rxn.products_[g].decay_rate_; score += simulation::keff * bank.wgt * rate * flux; @@ -879,8 +880,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -904,8 +905,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, continue; } else { score = 0.; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -943,27 +944,27 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->absorb_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_absorb_ * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->last_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_last_ * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } } else { @@ -975,8 +976,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * atom_density * flux; } } else { - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1002,8 +1003,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case ELASTIC: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Check if event MT matches - if (p->event_MT != ELASTIC) continue; - score = p->last_wgt * flux; + if (p->event_mt_ != ELASTIC) continue; + score = p->wgt_last_ * flux; } else { if (i_nuclide >= 0) { if (simulation::micro_xs[i_nuclide].elastic == CACHE_INVALID) @@ -1011,8 +1012,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = simulation::micro_xs[i_nuclide].elastic * atom_density * flux; } else { score = 0.; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1036,39 +1037,39 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->E_last_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score = p->absorb_wgt * q_value - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_absorb_ * q_value + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - const auto& nuc {*data::nuclides[p->event_nuclide]}; - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->E_last_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->E_last_); } - score = p->last_wgt * q_value - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->wgt_last_ * q_value + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } } else { @@ -1077,16 +1078,16 @@ score_general_ce(const Particle* p, int i_tally, int start_index, double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->E_last_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->E_last_); } score = q_value * simulation::micro_xs[i_nuclide].fission * atom_density * flux; } else { - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1094,10 +1095,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->E_last_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->E_last_); } score += q_value * simulation::micro_xs[j_nuclide].fission * atom_density * flux; @@ -1116,8 +1117,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case N_A: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Check if the event MT matches - if (p->event_MT != score_bin) continue; - score = p->last_wgt * flux; + if (p->event_mt_ != score_bin) continue; + score = p->wgt_last_ * flux; } else { int m; switch (score_bin) { @@ -1133,8 +1134,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * flux; } else { score = 0.; - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1151,8 +1152,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (tally.estimator_ == ESTIMATOR_ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need // to check if it matches the MT number of the event - if (p->event_MT != score_bin) continue; - score = p->last_wgt*flux; + if (p->event_mt_ != score_bin) continue; + score = p->wgt_last_*flux; } else { // Any other cross section has to be calculated on-the-fly if (score_bin < 2) fatal_error("Invalid score type on tally " @@ -1174,8 +1175,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } } else { - if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material]}; + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1223,7 +1224,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, auto& tally {*model::tallies[i_tally]}; // Set the direction and group to use with get_xs - const double* p_uvw; + Direction p_u; int p_g; if (tally.estimator_ == ESTIMATOR_ANALOG || tally.estimator_ == ESTIMATOR_COLLISION) { @@ -1232,41 +1233,41 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // Then we either are alive and had a scatter (and so g changed), // or are dead and g did not change - if (p->alive) { - p_uvw = p->last_uvw; - p_g = p->last_g; + if (p->alive_) { + p_u = p->u_last_; + p_g = p->g_last_; } else { - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_u = p->u_local(); + p_g = p->g_; } - } else if (p->event == EVENT_SCATTER) { + } else if (p->event_ == EVENT_SCATTER) { // Then the energy group has been changed by the scattering routine // meaning gin is now in p % last_g - p_uvw = p->last_uvw; - p_g = p->last_g; + p_u = p->u_last_; + p_g = p->g_last_; } else { // No scatter, no change in g. - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_u = p->u_local(); + p_g = p->g_; } } else { // No actual collision so g has not changed. - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_u = p->u_local(); + p_g = p->g_; } // To significantly reduce de-referencing, point matxs to the macroscopic // Mgxs for the material of interest - data::macro_xs[p->material].set_angle_index(p_uvw); + data::macro_xs[p->material_].set_angle_index(p_u); // Do same for nucxs, point it to the microscopic nuclide data of interest if (i_nuclide >= 0) { // And since we haven't calculated this temperature index yet, do so now - data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT); - data::nuclides_MG[i_nuclide].set_angle_index(p_uvw); + data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT_); + data::nuclides_MG[i_nuclide].set_angle_index(p_u); } for (auto i = 0; i < tally.scores_.size(); ++i) { @@ -1286,9 +1287,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->wgt_last_ + p->wgt_absorb_; } else { - score = p->last_wgt; + score = p->wgt_last_; } score *= flux / simulation::material_xs.total; } else { @@ -1304,15 +1305,15 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->wgt_last_ + p->wgt_absorb_; } else { - score = p->last_wgt; + score = p->wgt_last_; } //TODO: should flux be multiplied in above instead of below? if (i_nuclide >= 0) { score *= flux * atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } } else { if (i_nuclide >= 0) { @@ -1334,18 +1335,18 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->wgt_last_ + p->wgt_absorb_; } else { - score = p->last_wgt; + score = p->wgt_last_; } if (i_nuclide >= 0) { score *= flux * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } else { score *= flux - * get_macro_xs(p->material, MG_GET_XS_INVERSE_VELOCITY, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + * get_macro_xs(p->material_, MG_GET_XS_INVERSE_VELOCITY, p_g) + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } } else { if (i_nuclide >= 0) { @@ -1353,7 +1354,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_INVERSE_VELOCITY, p_g); + * get_macro_xs(p->material_, MG_GET_XS_INVERSE_VELOCITY, p_g); } } break; @@ -1362,24 +1363,24 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_SCATTER) continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, - p->last_g, &p->g, &p->mu, nullptr) - / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU_MULT, - p->last_g, &p->g, &p->mu, nullptr); + p->g_last_, &p->g_, &p->mu_, nullptr) + / get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU_MULT, + p->g_last_, &p->g_, &p->mu_, nullptr); } } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( - i_nuclide, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); + i_nuclide, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); } else { score = flux * get_macro_xs( - p->material, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); + p->material_, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); } } break; @@ -1388,27 +1389,27 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_NU_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_SCATTER) continue; // For scattering production, we need to use the pre-collision weight // times the multiplicity as the estimate for the number of neutrons // exiting a reaction with neutrons in the exit channel - score = p->wgt * flux; + score = p->wgt_ * flux; // Since we transport based on material data, the angle selected // was not selected from the f(mu) for the nuclide. Therefore // adjust the score by the actual probability for that nuclide. if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU, - p->last_g, &p->g, &p->mu, nullptr) - / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU, - p->last_g, &p->g, &p->mu, nullptr); + p->g_last_, &p->g_, &p->mu_, nullptr) + / get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU, + p->g_last_, &p->g_, &p->mu_, nullptr); } } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( i_nuclide, MG_GET_XS_SCATTER, p_g); } else { - score = flux * get_macro_xs(p->material, MG_GET_XS_SCATTER, p_g); + score = flux * get_macro_xs(p->material_, MG_GET_XS_SCATTER, p_g); } } break; @@ -1419,18 +1420,18 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { @@ -1449,30 +1450,30 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { score = get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_FISSION, p_g) * flux; + score = get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g) * flux; } } break; @@ -1480,7 +1481,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_NU_FISSION: if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1492,29 +1493,29 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p->wgt_bank * flux; + score = simulation::keff * p->wgt_bank_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } else { @@ -1522,7 +1523,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_NU_FISSION, p_g) * flux; + score = get_macro_xs(p->material_, MG_GET_XS_NU_FISSION, p_g) * flux; } } break; @@ -1530,7 +1531,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_PROMPT_NU_FISSION: if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1542,32 +1543,32 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_PROMPT_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p->n_delayed_bank, - p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank); - score = simulation::keff * p->wgt_bank * prompt_frac * flux; + auto n_delayed = std::accumulate(p->n_delayed_bank_, + p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank_); + score = simulation::keff * p->wgt_bank_ * prompt_frac * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } else { @@ -1575,7 +1576,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) + score = get_macro_xs(p->material_, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux; } } @@ -1584,7 +1585,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing || p->fission) { + if (settings::survival_biasing || p->fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1596,7 +1597,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g) > 0) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -1605,17 +1606,17 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1625,21 +1626,21 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1655,26 +1656,26 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - score = simulation::keff * p->wgt_bank / p->n_bank - * p->n_delayed_bank[d-1] * flux; + score = simulation::keff * p->wgt_bank_ / p->n_bank_ + * p->n_delayed_bank_[d-1] * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p->n_delayed_bank, - p->n_delayed_bank+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p->wgt_bank / p->n_bank * n_delayed + auto n_delayed = std::accumulate(p->n_delayed_bank_, + p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + score = simulation::keff * p->wgt_bank_ / p->n_bank_ * n_delayed * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } @@ -1693,7 +1694,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1705,7 +1706,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g); + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g); } } } @@ -1718,7 +1719,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g) > 0) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -1727,21 +1728,21 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1755,26 +1756,26 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = 0.; for (auto d = 0; d < data::num_delayed_groups; ++d) { if (i_nuclide >= 0) { - score += p->absorb_wgt * flux + score += p->wgt_absorb_ * flux * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { - score += p->absorb_wgt * flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + score += p->wgt_absorb_ * flux + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } } } } else { // Skip any non-fission events - if (!p->fission) continue; + if (!p->fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1783,8 +1784,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p->n_bank; ++i) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { @@ -1793,10 +1794,10 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &g) * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } else { score += simulation::keff * bank.wgt * flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, p_g, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &g); } if (tally.delayedgroup_filter_ != C_NONE) { @@ -1834,9 +1835,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score += flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1853,9 +1854,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score += flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } } @@ -1870,30 +1871,30 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - score = p->absorb_wgt * flux; + score = p->wgt_absorb_ * flux; } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->last_wgt * flux; + score = p->wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_KAPPA_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { score = get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) + score = get_macro_xs(p->material_, MG_GET_XS_KAPPA_FISSION, p_g) * flux; } } @@ -1923,7 +1924,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, int filter_index) { const Tally& tally {*model::tallies[i_tally]}; - const Material& material {*model::materials[p->material]}; + const Material& material {*model::materials[p->material_]}; // Score all individual nuclide reaction rates. for (auto i = 0; i < material.nuclide_.size(); ++i) { @@ -1980,7 +1981,7 @@ void score_analog_tally_ce(const Particle* p) // the event nuclide or the total material. Note that the i_nuclide // and flux arguments for score_general are not used for analog // tallies. - if (i_nuclide == p->event_nuclide || i_nuclide == -1) + if (i_nuclide == p->event_nuclide_ || i_nuclide == -1) score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); } @@ -1989,7 +1990,7 @@ void score_analog_tally_ce(const Particle* p) // In the case that the user has requested to tally all nuclides, we // can take advantage of the fact that we know exactly how nuclide // bins correspond to nuclide indices. First, tally the nuclide. - auto i = p->event_nuclide; + auto i = p->event_nuclide_; score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); @@ -2035,9 +2036,9 @@ void score_analog_tally_mg(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material]->atom_density_(j); + atom_density = model::materials[p->material_]->atom_density_(j); } score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, @@ -2061,7 +2062,7 @@ void score_tracklength_tally(const Particle* p, double distance) { // Determine the tracklength estimate of the flux - double flux = p->wgt * distance; + double flux = p->wgt_ * distance; for (auto i_tally : model::active_tracklength_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2080,7 +2081,7 @@ score_tracklength_tally(const Particle* p, double distance) // Loop over nuclide bins. if (tally.all_nuclides_) { - if (p->material != MATERIAL_VOID) + if (p->material_ != MATERIAL_VOID) score_all_nuclides(p, i_tally, flux*filter_weight, filter_index); } else { @@ -2089,10 +2090,10 @@ score_tracklength_tally(const Particle* p, double distance) double atom_density = 0.; if (i_nuclide >= 0) { - if (p->material != MATERIAL_VOID) { - auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; + if (p->material_ != MATERIAL_VOID) { + auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material]->atom_density_(j); + atom_density = model::materials[p->material_]->atom_density_(j); } } @@ -2126,9 +2127,9 @@ void score_collision_tally(const Particle* p) // Determine the collision estimate of the flux double flux; if (!settings::survival_biasing) { - flux = p->last_wgt / simulation::material_xs.total; + flux = p->wgt_last_ / simulation::material_xs.total; } else { - flux = (p->last_wgt + p->absorb_wgt) / simulation::material_xs.total; + flux = (p->wgt_last_ + p->wgt_absorb_) / simulation::material_xs.total; } for (auto i_tally : model::active_collision_tallies) { @@ -2156,9 +2157,9 @@ void score_collision_tally(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material]->atom_density_(j); + atom_density = model::materials[p->material_]->atom_density_(j); } //TODO: consider replacing this "if" with pointers or templates @@ -2189,7 +2190,7 @@ void score_surface_tally(const Particle* p, const std::vector& tallies) { // No collision, so no weight change when survival biasing - double flux = p->wgt; + double flux = p->wgt_; for (auto i_tally : tallies) { auto& tally {*model::tallies[i_tally]}; diff --git a/src/track_output.cpp b/src/track_output.cpp index df5cef2fc..50d6d1ef6 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -40,14 +40,14 @@ void add_particle_track() void write_particle_track(const Particle& p) { - tracks.back().push_back({p.coord[0].xyz}); + tracks.back().push_back(p.r()); } void finalize_particle_track(const Particle& p) { std::stringstream filename; filename << settings::path_output << "track_" << simulation::current_batch - << '_' << simulation::current_gen << '_' << p.id << ".h5"; + << '_' << simulation::current_gen << '_' << p.id_ << ".h5"; // Determine number of coordinates for each particle std::vector n_coords; diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 0f18e245f..93b729bd6 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -96,7 +96,6 @@ std::vector VolumeCalculation::execute() const std::vector> indices(n); std::vector> hits(n); Particle p; - p.initialize(); prn_set_stream(STREAM_VOLUME); @@ -105,41 +104,37 @@ std::vector VolumeCalculation::execute() const for (int i = i_start; i < i_end; i++) { set_particle_seed(i); - p.n_coord = 1; + p.n_coord_ = 1; Position xi {prn(), prn(), prn()}; - Position r {lower_left_ + xi*(upper_right_ - lower_left_)}; - // TODO: assign directly when xyz is Position - std::copy(&r.x, &r.x + 3, p.coord[0].xyz); - p.coord[0].uvw[0] = 0.5; - p.coord[1].uvw[1] = 0.5; - p.coord[2].uvw[2] = 0.5; + p.r() = lower_left_ + xi*(upper_right_ - lower_left_); + p.u() = {0.5, 0.5, 0.5}; // If this location is not in the geometry at all, move on to next block if (!find_cell(&p, false)) continue; if (domain_type_ == FILTER_MATERIAL) { - if (p.material != MATERIAL_VOID) { + if (p.material_ != MATERIAL_VOID) { for (int i_domain = 0; i_domain < n; i_domain++) { - if (model::materials[p.material]->id_ == domain_ids_[i_domain]) { - this->check_hit(p.material, indices[i_domain], hits[i_domain]); + if (model::materials[p.material_]->id_ == domain_ids_[i_domain]) { + this->check_hit(p.material_, indices[i_domain], hits[i_domain]); break; } } } } else if (domain_type_ == FILTER_CELL) { - for (int level = 0; level < p.n_coord; ++level) { + for (int level = 0; level < p.n_coord_; ++level) { for (int i_domain=0; i_domain < n; i_domain++) { - if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) { - this->check_hit(p.material, indices[i_domain], hits[i_domain]); + if (model::cells[p.coord_[level].cell]->id_ == domain_ids_[i_domain]) { + this->check_hit(p.material_, indices[i_domain], hits[i_domain]); break; } } } } else if (domain_type_ == FILTER_UNIVERSE) { - for (int level = 0; level < p.n_coord; ++level) { + for (int level = 0; level < p.n_coord_; ++level) { for (int i_domain = 0; i_domain < n; ++i_domain) { - if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) { - check_hit(p.material, indices[i_domain], hits[i_domain]); + if (model::universes[p.coord_[level].universe]->id_ == domain_ids_[i_domain]) { + check_hit(p.material_, indices[i_domain], hits[i_domain]); break; } } diff --git a/tests/regression_tests/sourcepoint_batch/test.py b/tests/regression_tests/sourcepoint_batch/test.py index 3086b5610..fdc5ecca1 100644 --- a/tests/regression_tests/sourcepoint_batch/test.py +++ b/tests/regression_tests/sourcepoint_batch/test.py @@ -19,7 +19,7 @@ class SourcepointTestHarness(TestHarness): # Read the statepoint file. with StatePoint(self._sp_name) as sp: # Add the source information. - xyz = sp.source[0]['xyz'] + xyz = sp.source[0]['r'] outstr += ' '.join(['{0:12.6E}'.format(x) for x in xyz]) outstr += "\n" diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index d5943b7be..0ebc898b4 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -236,6 +236,7 @@ def test_source_bank(capi_run): source = openmc.capi.source_bank() assert np.all(source['E'] > 0.0) assert np.all(source['wgt'] == 1.0) + assert np.allclose(np.linalg.norm(source['u'], axis=1), 1.0) def test_by_batch(capi_run):