From 368f89697dcd0fb19f41caaf5462644377f40131 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Mar 2019 14:46:21 -0600 Subject: [PATCH] Use Position and Direction in Particle class --- include/openmc/bank.h | 22 +-- include/openmc/capi.h | 18 +-- include/openmc/hdf5_interface.h | 10 ++ include/openmc/math_functions.h | 4 +- include/openmc/mesh.h | 2 +- include/openmc/mgxs.h | 8 +- include/openmc/mgxs_interface.h | 2 +- include/openmc/particle.h | 51 +++++-- include/openmc/physics.h | 14 +- include/openmc/physics_mg.h | 2 +- include/openmc/position.h | 5 +- include/openmc/source.h | 5 +- openmc/capi/core.py | 6 +- openmc/cmfd.py | 8 +- src/bank.cpp | 12 +- src/bremsstrahlung.cpp | 2 +- src/eigenvalue.cpp | 4 +- src/geometry.cpp | 89 +++++------- src/initialize.cpp | 12 +- src/math_functions.cpp | 68 +++++---- src/mesh.cpp | 16 +-- src/mgxs.cpp | 22 +-- src/mgxs_interface.cpp | 4 +- src/output.cpp | 8 +- src/particle.cpp | 101 +++++--------- src/particle_restart.cpp | 13 +- src/photon.cpp | 22 +-- src/physics.cpp | 132 ++++++++---------- src/physics_mg.cpp | 18 ++- src/plot.cpp | 118 +++++++--------- src/position.cpp | 6 + src/source.cpp | 22 ++- src/state_point.cpp | 24 ++-- src/tallies/derivative.cpp | 4 +- src/tallies/filter_azimuthal.cpp | 4 +- src/tallies/filter_mesh.cpp | 2 +- src/tallies/filter_polar.cpp | 4 +- src/tallies/filter_sph_harm.cpp | 2 +- src/tallies/filter_sptl_legendre.cpp | 6 +- src/tallies/filter_zernike.cpp | 8 +- src/tallies/tally_scoring.cpp | 17 +-- src/track_output.cpp | 2 +- src/volume_calc.cpp | 8 +- .../sourcepoint_batch/test.py | 2 +- tests/unit_tests/test_capi.py | 1 + 45 files changed, 423 insertions(+), 487 deletions(-) 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/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/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 f753cb119..0c3e8738f 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,15 +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}; +//============================================================================== +// 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 @@ -59,6 +63,16 @@ public: neutron, photon, electron, positron }; + //! Saved ("banked") state of a particle + struct Bank { + Position r; + Direction u; + double E; + double wgt; + int delayed_group; + Type particle; + }; + // Constructors Particle(); @@ -85,13 +99,13 @@ public: 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_ {1.0}; //!< pre-collision particle weight - double absorb_wgt_ {0.0}; //!< 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 last_wgt_ {1.0}; //!< pre-collision particle weight + double absorb_wgt_ {0.0}; //!< weight absorbed for survival biasing // What event took place bool fission_ {false}; //!< did particle cause implicit fission @@ -126,6 +140,18 @@ public: int64_t n_secondary_ {}; Bank secondary_bank_[MAX_SECONDARY]; + Position& r() { return coord_[0].r; } + const Position& r() const { return coord_[0].r; } + + Position& r_local() { return coord_[n_coord_ - 1].r; } + const Position& r_local() const { return coord_[n_coord_ - 1].r; } + + Direction& u() { return coord_[0].u; } + const Direction& u() const { return coord_[0].u; } + + 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(); @@ -133,11 +159,10 @@ public: // //! 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, Type type, bool run_CE); + 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/source.h b/include/openmc/source.h index bac871d6d..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,7 +39,7 @@ public: //! Sample from the external source distribution //! \return Sampled site - Bank sample() const; + Particle::Bank sample() const; // Properties double strength() const { return strength_; } @@ -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/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 931b95faa..3b9d8b591 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -108,7 +108,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon - p.create_secondary(p.coord_[0].uvw, w, Particle::Type::photon, true); + p.create_secondary(p.u(), w, Particle::Type::photon); *E_lost += w; } } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 1cf159889..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 f9aea182e..4b0dc9184 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -44,7 +44,7 @@ check_cell_overlap(Particle* p) // 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 (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_ << ", " @@ -78,8 +78,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) 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}; + 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; @@ -99,8 +99,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) 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}; + 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; @@ -169,36 +169,28 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) 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]; + 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; } @@ -213,18 +205,12 @@ 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}; - auto i_xyz = lat.get_indices(r, u); + 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_; @@ -324,10 +310,7 @@ cross_lattice(Particle* p, int lattice_translation[3]) 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. @@ -377,8 +360,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, // 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}; + 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. @@ -399,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; } @@ -466,15 +449,13 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) { Particle p; - 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; } 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/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 4f0a746c7..a72fce356 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; @@ -522,9 +522,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_; @@ -652,7 +652,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 @@ -670,7 +670,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) { 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 9bf785a7a..c30850e8f 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -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/output.cpp b/src/output.cpp index 8d41a9889..92ad14e9d 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -182,10 +182,10 @@ extern "C" void print_particle(Particle* p) << 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. diff --git a/src/particle.cpp b/src/particle.cpp index 059d4ac41..d5831b1d4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -66,19 +66,18 @@ Particle::clear() } void -Particle::create_secondary(const double* uvw, double E, Type type, bool run_CE) +Particle::create_secondary(Direction u, double E, Type type) { if (n_secondary_ == MAX_SECONDARY) { fatal_error("Too many secondary particles created."); } int64_t n = n_secondary_; - secondary_bank_[n].particle = static_cast(type); + secondary_bank_[n].particle = type; secondary_bank_[n].wgt = wgt_; - std::copy(coord_[0].xyz, coord_[0].xyz + 3, secondary_bank_[n].xyz); - std::copy(uvw, uvw + 3, secondary_bank_[n].uvw); - secondary_bank_[n].E = E; - if (!run_CE) secondary_bank_[n].E = g_; + secondary_bank_[n].r = this->r(); + secondary_bank_[n].u = u; + secondary_bank_[n].E = settings::run_CE ? E : g_; ++n_secondary_; } @@ -95,14 +94,14 @@ Particle::from_source(const Bank* src) fission_ = false; // copy attributes from source bank site - type_ = static_cast(src->particle); + 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_); + 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; @@ -153,8 +152,8 @@ Particle::transport() // 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_); + 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 @@ -186,7 +185,7 @@ Particle::transport() } } 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); @@ -225,10 +224,7 @@ Particle::transport() // 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]; + coord_[j].r += distance * coord_[j].u; } // Score track-length tallies @@ -323,25 +319,25 @@ Particle::transport() 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; - // Set all uvws to base level -- right now, after a collision, only the - // base level uvws are changed + // 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& 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; } } @@ -403,10 +399,7 @@ 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); } @@ -443,22 +436,17 @@ 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]; @@ -476,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) { @@ -502,12 +488,10 @@ 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 @@ -516,11 +500,7 @@ Particle::cross_surface() 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 @@ -538,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) { @@ -588,9 +566,7 @@ Particle::cross_surface() // 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]; + this->r() += TINY_BIT * this->u(); // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. @@ -673,9 +649,8 @@ Particle::write_restart() const 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 7c0772879..a188d1104 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -45,11 +45,8 @@ void read_particle_restart(Particle& p, int& previous_run_mode) p.type_ = static_cast(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, "xyz", p.r()); + read_dataset(file_id, "uvw", p.u()); // Set energy group and average energy in multi-group mode if (!settings::run_CE) { @@ -59,9 +56,9 @@ void read_particle_restart(Particle& p, int& previous_run_mode) // 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.r_last_current_ = p.r(); + p.r_last_ = p.r(); + p.u_last_ = p.u(); p.last_E_ = p.E_; p.last_g_ = p.g_; diff --git a/src/photon.cpp b/src/photon.cpp index 44ff39508..2dd0928ba 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -653,12 +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; - p.create_secondary(uvw.data(), E, Particle::Type::photon, true); + p.create_secondary(u, E, Particle::Type::photon); return; } @@ -678,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); @@ -690,7 +690,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-radiative transition -- Auger/Coster-Kronig effect // Create auger electron - p.create_secondary(uvw.data(), E, Particle::Type::electron, true); + p.create_secondary(u, E, Particle::Type::electron); // Fill hole left by emitted auger electron int i_hole = shell_map_.at(secondary); @@ -700,7 +700,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Radiative transition -- get X-ray energy // Create fluorescent photon - p.create_secondary(uvw.data(), E, Particle::Type::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 28a795254..39f73b235 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -136,8 +136,8 @@ void sample_neutron_reaction(Particle* p) } 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 @@ -181,12 +181,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ 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(Particle::Type::neutron); + bank_array[i].particle = Particle::Type::neutron; // Set the weight of the fission bank site bank_array[i].wgt = 1. / weight; @@ -244,7 +242,7 @@ 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->u() = rotate_angle(p->u(), mu, nullptr); p->event_mt_ = COHERENT; return; } @@ -270,10 +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); - p->create_secondary(uvw, E_electron, Particle::Type::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 @@ -285,7 +281,7 @@ 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->u() = rotate_angle(p->u(), mu, &phi); p->event_mt_ = INCOHERENT; return; } @@ -326,13 +322,13 @@ 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 - p->create_secondary(uvw.data(), E_electron, Particle::Type::electron, true); + p->create_secondary(u, E_electron, Particle::Type::electron); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons @@ -355,15 +351,12 @@ 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); - p->create_secondary(uvw, E_electron, Particle::Type::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); - p->create_secondary(uvw, E_positron, Particle::Type::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; @@ -396,18 +389,14 @@ 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 - p->create_secondary(uvw.data(), MASS_ELECTRON_EV, Particle::Type::photon, true); - - uvw[0] = -uvw[0]; - uvw[1] = -uvw[1]; - uvw[2] = -uvw[2]; - p->create_secondary(uvw.data(), MASS_ELECTRON_EV, Particle::Type::photon, true); + p->create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon); + p->create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon); p->E_ = 0.0; p->alive_ = false; @@ -585,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]}; @@ -613,8 +602,8 @@ void scatter(Particle* p, int i_nuclide) 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; sampled = true; @@ -625,7 +614,7 @@ 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; sampled = true; @@ -673,39 +662,33 @@ void scatter(Particle* p, int i_nuclide) // 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 @@ -720,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; } @@ -739,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]}; @@ -766,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; @@ -973,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 @@ -983,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]}; @@ -1098,14 +1078,14 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) 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) { - p->create_secondary(p->coord_[0].uvw, p->E_, Particle::Type::neutron, true); + p->create_secondary(p->u(), p->E_, Particle::Type::neutron); } } else { // Otherwise, change weight of particle based on yield @@ -1135,12 +1115,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide) 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 - p->create_secondary(uvw, E, Particle::Type::photon, true); + p->create_secondary(u, E, Particle::Type::photon); } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index bfccbb589..b9d1d8279 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -92,7 +92,7 @@ scatter(Particle* p) 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]; @@ -102,7 +102,7 @@ scatter(Particle* p) } 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 @@ -152,12 +152,10 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, 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(Particle::Type::neutron); + bank_array[i].particle = Particle::Type::neutron; // Set the weight of the fission bank site bank_array[i].wgt = 1. / weight; @@ -168,16 +166,16 @@ 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); + bank_array[i].E = gout + 1; bank_array[i].delayed_group = dg + 1; // Set the delayed group on the particle as well diff --git a/src/plot.cpp b/src/plot.cpp index c3456b3ae..ae2f00023 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -106,48 +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; - std::copy(xyz, xyz+3, p.coord_[0].xyz); - std::copy(dir, dir+3, p.coord_[0].uvw); + 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; } @@ -308,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 " @@ -328,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 " @@ -339,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 " @@ -745,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++) { @@ -782,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 @@ -845,16 +832,13 @@ 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; - std::copy(ll.begin(), ll.begin()+ll.size(), p.coord_[0].xyz); - std::copy(dir, dir+3, p.coord_[0].uvw); + p.r() = ll; + p.u() = u; p.coord_[0].universe = model::root_universe; // Open binary plot file for writing @@ -889,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]]; @@ -908,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/source.cpp b/src/source.cpp index ed1c89b8d..6b24270a4 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -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,14 @@ 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(); // 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(&site.r.x, &cell_index, &instance); found = (err != OPENMC_E_GEOMETRY); // Check if spatial site is in fissionable material @@ -200,10 +197,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 +284,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 +306,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/tallies/derivative.cpp b/src/tallies/derivative.cpp index 1ddd7af83..b9a1a82a1 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -667,8 +667,8 @@ void score_collision_derivative(const Particle* p) 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 4b78cabd5..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_mesh.cpp b/src/tallies/filter_mesh.cpp index 518470bf4..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_polar.cpp b/src/tallies/filter_polar.cpp index 20eb7608b..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 3ef2df43c..9aae06bee 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -44,7 +44,7 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator, // 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 14b8fe14e..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_zernike.cpp b/src/tallies/filter_zernike.cpp index 93ff16109..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_scoring.cpp b/src/tallies/tally_scoring.cpp index 1931c37ce..6b3769222 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" @@ -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) { @@ -1233,40 +1234,40 @@ 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_u = p->u_last_; p_g = p->last_g_; } else { - p_uvw = p->coord_[p->n_coord_-1].uvw; + p_u = p->u_local(); p_g = p->g_; } } 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_u = p->u_last_; p_g = p->last_g_; } else { // No scatter, no change in g. - p_uvw = p->coord_[p->n_coord_-1].uvw; + 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_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_angle_index(p_u); } for (auto i = 0; i < tally.scores_.size(); ++i) { diff --git a/src/track_output.cpp b/src/track_output.cpp index 34cf7fdf6..50d6d1ef6 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -40,7 +40,7 @@ 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) diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index a23481539..93b729bd6 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -106,12 +106,8 @@ std::vector VolumeCalculation::execute() const 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; 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):