crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Geometron (#2744)

Browse source
This commit is contained in:
Gavin Ridley 2024-01-16 10:35:57 -06:00 • committed by GitHub
parent 971c5f77a5
commit 5549b58e1f
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
16 changed files with 403 additions and 284 deletions
Download patch file Download diff file Expand all files Collapse all files

14
include/openmc/cell.h
View file

@ -40,6 +40,7 @@ constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
//==============================================================================
class Cell;
class GeometryState;
class ParentCell;
class CellInstance;
class Universe;
@ -82,7 +83,7 @@ public:
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const;
Position r, Direction u, int32_t on_surface) const;
//! Get the BoundingBox for this cell.
BoundingBox bounding_box(int32_t cell_id) const;
@ -183,7 +184,7 @@ public:
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const = 0;
Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
@ -260,7 +261,8 @@ protected:
//! \param[in] instance of the cell to find parent cells for
//! \param[in] p particle used to do a fast search for parent cells
//! \return parent cells
vector<ParentCell> find_parent_cells(int32_t instance, Particle& p) const;
vector<ParentCell> find_parent_cells(
int32_t instance, GeometryState& p) const;
//! Determine the path to this cell instance in the geometry hierarchy
//! \param[in] instance of the cell to find parent cells for
@ -332,10 +334,10 @@ public:
// Methods
vector<int32_t> surfaces() const override { return region_.surfaces(); }
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const override
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override
{
return region_.distance(r, u, on_surface, p);
return region_.distance(r, u, on_surface);
}
bool contains(Position r, Direction u, int32_t on_surface) const override

8
include/openmc/dagmc.h
View file

@ -42,7 +42,7 @@ public:
double evaluate(Position r) const override;
double distance(Position r, Direction u, bool coincident) const override;
Direction normal(Position r) const override;
Direction reflect(Position r, Direction u, Particle* p) const override;
Direction reflect(Position r, Direction u, GeometryState* p) const override;
inline void to_hdf5_inner(hid_t group_id) const override {};
@ -63,8 +63,8 @@ public:
bool contains(Position r, Direction u, int32_t on_surface) const override;
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const override;
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override;
BoundingBox bounding_box() const override;
@ -143,7 +143,7 @@ public:
//! string of the ID ranges for entities of dimension \p dim
std::string dagmc_ids_for_dim(int dim) const;
bool find_cell(Particle& p) const override;
bool find_cell(GeometryState& p) const override;
void to_hdf5(hid_t universes_group) const override;

16
include/openmc/geometry.h
View file

@ -11,7 +11,7 @@
namespace openmc {
class BoundaryInfo;
class Particle;
class GeometryState;
//==============================================================================
// Global variables
@ -39,7 +39,7 @@ inline bool coincident(double d1, double d2)
//! Check for overlapping cells at a particle's position.
//==============================================================================
bool check_cell_overlap(Particle& p, bool error = true);
bool check_cell_overlap(GeometryState& p, bool error = true);
//==============================================================================
//! Get the cell instance for a particle at the specified universe level
@ -50,7 +50,7 @@ bool check_cell_overlap(Particle& p, bool error = true);
//! should be computed. \return The instance of the cell at the specified level.
//==============================================================================
int cell_instance_at_level(const Particle& p, int level);
int cell_instance_at_level(const GeometryState& p, int level);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -60,20 +60,22 @@ int cell_instance_at_level(const Particle& p, int level);
//! \return True if the particle's location could be found and ascribed to a
//! valid geometry coordinate stack.
//==============================================================================
bool exhaustive_find_cell(Particle& p);
bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings
bool exhaustive_find_cell(GeometryState& p, bool verbose = false);
bool neighbor_list_find_cell(
GeometryState& p, bool verbose = false); // Only usable on surface crossings
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
void cross_lattice(Particle& p, const BoundaryInfo& boundary);
void cross_lattice(
GeometryState& p, const BoundaryInfo& boundary, bool verbose = false);
//==============================================================================
//! Find the next boundary a particle will intersect.
//==============================================================================
BoundaryInfo distance_to_boundary(Particle& p);
BoundaryInfo distance_to_boundary(GeometryState& p);
} // namespace openmc

14
include/openmc/particle.h
View file

@ -5,7 +5,6 @@
//! \brief Particle type
#include <cstdint>
#include <sstream>
#include <string>
#include "openmc/constants.h"
@ -103,17 +102,8 @@ public:
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message)
{
mark_as_lost(message.c_str());
}
void mark_as_lost(const std::stringstream& message)
{
mark_as_lost(message.str());
}
virtual void mark_as_lost(const char* message) override;
using GeometryState::mark_as_lost;
//! create a particle restart HDF5 file
void write_restart() const;

440
include/openmc/particle_data.h
View file

@ -194,6 +194,157 @@ struct BoundaryInfo {
lattice_translation {}; //!< which way lattice indices will change
};
/*
* Contains all geometry state information for a particle.
*/
class GeometryState {
public:
GeometryState();
/*
* GeometryState does not store any ID info, so give some reasonable behavior
* here. The Particle class redefines this. This is only here for the error
* reporting behavior that occurs in geometry.cpp. The explanation for
* mark_as_lost is the same.
*/
virtual void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message);
void mark_as_lost(const std::stringstream& message);
// resets all coordinate levels for the particle
void clear()
{
for (auto& level : coord_)
level.reset();
n_coord_ = 1;
}
// Initialize all internal state from position and direction
void init_from_r_u(Position r_a, Direction u_a)
{
clear();
surface() = 0;
material() = C_NONE;
r() = r_a;
u() = u_a;
r_last_current() = r_a;
r_last() = r_a;
u_last() = u_a;
}
// Unique ID. This is not geometric info, but the
// error reporting in geometry.cpp requires this.
// We could save this to implement it in Particle,
// but that would require virtuals.
int64_t& id() { return id_; }
const int64_t& id() const { return id_; }
// Number of current coordinate levels
int& n_coord() { return n_coord_; }
const int& n_coord() const { return n_coord_; }
// Offset for distributed properties
int& cell_instance() { return cell_instance_; }
const int& cell_instance() const { return cell_instance_; }
// Coordinates for all nesting levels
LocalCoord& coord(int i) { return coord_[i]; }
const LocalCoord& coord(int i) const { return coord_[i]; }
const vector<LocalCoord>& coord() const { return coord_; }
// Innermost universe nesting coordinates
LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; }
const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; }
// Last coordinates on all nesting levels, before crossing a surface
int& n_coord_last() { return n_coord_last_; }
const int& n_coord_last() const { return n_coord_last_; }
int& cell_last(int i) { return cell_last_[i]; }
const int& cell_last(int i) const { return cell_last_[i]; }
// Coordinates of last collision or reflective/periodic surface
// crossing for current tallies
Position& r_last_current() { return r_last_current_; }
const Position& r_last_current() const { return r_last_current_; }
// Previous direction and spatial coordinates before a collision
Position& r_last() { return r_last_; }
const Position& r_last() const { return r_last_; }
Position& u_last() { return u_last_; }
const Position& u_last() const { return u_last_; }
// 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; }
// Surface that the particle is on
int& surface() { return surface_; }
const int& surface() const { return surface_; }
// Boundary information
BoundaryInfo& boundary() { return boundary_; }
#ifdef DAGMC
// DagMC state variables
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
#endif
// material of current and last cell
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
const int& material_last() const { return material_last_; }
// temperature of current and last cell
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
private:
int64_t id_ {-1}; //!< Unique ID
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
vector<LocalCoord> coord_; //!< coordinates for all levels
int n_coord_last_ {1}; //!< number of current coordinates
vector<int> cell_last_; //!< coordinates for all levels
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
int surface_ {0}; //!< index for surface particle is on
BoundaryInfo boundary_; //!< Info about the next intersection
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
#ifdef DAGMC
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
};
//============================================================================
//! Defines how particle data is laid out in memory
//============================================================================
@ -229,163 +380,112 @@ struct BoundaryInfo {
* Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
* https://doi.org/10.1016/j.anucene.2017.11.032.
*/
class ParticleData {
public:
//----------------------------------------------------------------------------
// Constructors
ParticleData();
class ParticleData : public GeometryState {
private:
//==========================================================================
// Data members (accessor methods are below)
// Data members -- see public: below for descriptions
// Cross section caches
vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
vector<ElementMicroXS> photon_xs_; //!< Microscopic photon cross sections
MacroXS macro_xs_; //!< Macroscopic cross sections
CacheDataMG mg_xs_cache_; //!< Multigroup XS cache
vector<NuclideMicroXS> neutron_xs_;
vector<ElementMicroXS> photon_xs_;
MacroXS macro_xs_;
CacheDataMG mg_xs_cache_;
int64_t id_; //!< Unique ID
ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.)
ParticleType type_ {ParticleType::neutron};
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
vector<LocalCoord> coord_; //!< coordinates for all levels
double E_;
double E_last_;
int g_ {0};
int g_last_;
// Particle coordinates before crossing a surface
int n_coord_last_ {1}; //!< number of current coordinates
vector<int> cell_last_; //!< coordinates for all levels
double wgt_ {1.0};
double mu_;
double time_ {0.0};
double time_last_ {0.0};
double wgt_last_ {1.0};
// Energy data
double E_; //!< post-collision energy in eV
double E_last_; //!< pre-collision energy in eV
int g_ {C_NONE}; //!< post-collision energy group (MG only)
int g_last_; //!< pre-collision energy group (MG only)
bool fission_ {false};
TallyEvent event_;
int event_nuclide_;
int event_mt_;
int delayed_group_ {0};
// Other physical data
double wgt_ {1.0}; //!< particle weight
double mu_; //!< angle of scatter
double time_ {0.0}; //!< time in [s]
double time_last_ {0.0}; //!< previous time in [s]
int n_bank_ {0};
int n_bank_second_ {0};
double wgt_bank_ {0.0};
int n_delayed_bank_[MAX_DELAYED_GROUPS];
// Other physical data
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
int cell_born_ {-1};
// What event took place
bool fission_ {false}; //!< did particle cause implicit fission
TallyEvent event_; //!< scatter, absorption
int event_nuclide_; //!< index in nuclides array
int event_mt_; //!< reaction MT
int delayed_group_ {0}; //!< delayed group
int n_collision_ {0};
// Post-collision physical data
int n_bank_ {0}; //!< number of fission sites banked
int n_bank_second_ {0}; //!< number of secondary particles 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_ {0}; //!< index for surface particle is on
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
// Boundary information
BoundaryInfo boundary_;
// Temperature of current cell
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
// Statistical data
int n_collision_ {0}; //!< number of collisions
// Track output
bool write_track_ {false};
// Current PRNG state
uint64_t seeds_[N_STREAMS]; // current seeds
int stream_; // current RNG stream
uint64_t seeds_[N_STREAMS];
int stream_;
// Secondary particle bank
vector<SourceSite> secondary_bank_;
int64_t current_work_; // current work index
int64_t current_work_;
vector<double> flux_derivs_; // for derivatives for this particle
vector<double> flux_derivs_;
vector<FilterMatch> filter_matches_; // tally filter matches
vector<FilterMatch> filter_matches_;
vector<TrackStateHistory> tracks_; // tracks for outputting to file
vector<TrackStateHistory> tracks_;
vector<NuBank> nu_bank_; // bank of most recently fissioned particles
vector<NuBank> nu_bank_;
vector<double> pht_storage_; // interim pulse-height results
vector<double> pht_storage_;
// Global tally accumulators
double keff_tally_absorption_ {0.0};
double keff_tally_collision_ {0.0};
double keff_tally_tracklength_ {0.0};
double keff_tally_leakage_ {0.0};
bool trace_ {false}; //!< flag to show debug information
bool trace_ {false};
double collision_distance_; // distance to particle's next closest collision
double collision_distance_;
int n_event_ {0}; // number of events executed in this particle's history
int n_event_ {0};
// Weight window information
int n_split_ {0}; // Number of times this particle has been split
double ww_factor_ {
0.0}; // Particle-specific factor for on-the-fly weight window adjustment
int n_split_ {0};
double ww_factor_ {0.0};
// DagMC state variables
#ifdef DAGMC
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
int64_t n_progeny_ {0}; // Number of progeny produced by this particle
int64_t n_progeny_ {0};
public:
//----------------------------------------------------------------------------
// Constructors
ParticleData();
//==========================================================================
// Methods and accessors
NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; }
// Cross section caches
NuclideMicroXS& neutron_xs(int i)
{
return neutron_xs_[i];
} // Microscopic neutron cross sections
const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; }
// Microscopic photon cross sections
ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; }
// Macroscopic cross sections
MacroXS& macro_xs() { return macro_xs_; }
const MacroXS& macro_xs() const { return macro_xs_; }
// Multigroup macroscopic cross sections
CacheDataMG& mg_xs_cache() { return mg_xs_cache_; }
const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; }
int64_t& id() { return id_; }
const int64_t& id() const { return id_; }
// Particle type (n, p, e, gamma, etc)
ParticleType& type() { return type_; }
const ParticleType& type() const { return type_; }
int& n_coord() { return n_coord_; }
const int& n_coord() const { return n_coord_; }
int& cell_instance() { return cell_instance_; }
const int& cell_instance() const { return cell_instance_; }
LocalCoord& coord(int i) { return coord_[i]; }
const LocalCoord& coord(int i) const { return coord_[i]; }
const vector<LocalCoord>& coord() const { return coord_; }
LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; }
const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; }
int& n_coord_last() { return n_coord_last_; }
const int& n_coord_last() const { return n_coord_last_; }
int& cell_last(int i) { return cell_last_[i]; }
const int& cell_last(int i) const { return cell_last_[i]; }
// Current particle energy, energy before collision,
// and corresponding multigroup group indices. Energy
// units are eV.
double& E() { return E_; }
const double& E() const { return E_; }
double& E_last() { return E_last_; }
@ -395,113 +495,119 @@ public:
int& g_last() { return g_last_; }
const int& g_last() const { return g_last_; }
// Statistic weight of particle. Setting to zero
// indicates that the particle is dead.
double& wgt() { return wgt_; }
double wgt() const { return wgt_; }
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
bool alive() const { return wgt_ != 0.0; }
// Polar scattering angle after a collision
double& mu() { return mu_; }
const double& mu() const { return mu_; }
// Tracks the time of a particle as it traverses the problem.
// Units are seconds.
double& time() { return time_; }
const double& time() const { return time_; }
double& time_last() { return time_last_; }
const double& time_last() const { return time_last_; }
bool alive() const { return wgt_ != 0.0; }
Position& r_last_current() { return r_last_current_; }
const Position& r_last_current() const { return r_last_current_; }
Position& r_last() { return r_last_; }
const Position& r_last() const { return r_last_; }
Position& u_last() { return u_last_; }
const Position& u_last() const { return u_last_; }
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
bool& fission() { return fission_; }
// What event took place, described in greater detail below
TallyEvent& event() { return event_; }
const TallyEvent& event() const { return event_; }
int& event_nuclide() { return event_nuclide_; }
bool& fission() { return fission_; } // true if implicit fission
int& event_nuclide() { return event_nuclide_; } // index of collision nuclide
const int& event_nuclide() const { return event_nuclide_; }
int& event_mt() { return event_mt_; }
int& delayed_group() { return delayed_group_; }
int& event_mt() { return event_mt_; } // MT number of collision
int& delayed_group() { return delayed_group_; } // delayed group
int& n_bank() { return n_bank_; }
int& n_bank_second() { return n_bank_second_; }
double& wgt_bank() { return wgt_bank_; }
int* n_delayed_bank() { return n_delayed_bank_; }
int& n_delayed_bank(int i) { return n_delayed_bank_[i]; }
// Post-collision data
int& n_bank() { return n_bank_; } // number of banked fission sites
int& n_bank_second()
{
return n_bank_second_;
} // number of secondaries banked
double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites
int* n_delayed_bank()
{
return n_delayed_bank_;
} // number of delayed fission sites
int& n_delayed_bank(int i)
{
return n_delayed_bank_[i];
} // number of delayed fission sites
int& surface() { return surface_; }
const int& surface() const { return surface_; }
// Index of cell particle is born in
int& cell_born() { return cell_born_; }
const int& cell_born() const { return cell_born_; }
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
const int& material_last() const { return material_last_; }
BoundaryInfo& boundary() { return boundary_; }
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
// index of the current and last material
// Total number of collisions suffered by particle
int& n_collision() { return n_collision_; }
const int& n_collision() const { return n_collision_; }
// whether this track is to be written
bool& write_track() { return write_track_; }
// RNG state
uint64_t& seeds(int i) { return seeds_[i]; }
uint64_t* seeds() { return seeds_; }
int& stream() { return stream_; }
// secondary particle bank
SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
// Current simulation work index
int64_t& current_work() { return current_work_; }
const int64_t& current_work() const { return current_work_; }
// Used in tally derivatives
double& flux_derivs(int i) { return flux_derivs_[i]; }
const double& flux_derivs(int i) const { return flux_derivs_[i]; }
// Matches of tallies
decltype(filter_matches_)& filter_matches() { return filter_matches_; }
FilterMatch& filter_matches(int i) { return filter_matches_[i]; }
// Tracks to output to file
decltype(tracks_)& tracks() { return tracks_; }
// Bank of recently fissioned particles
decltype(nu_bank_)& nu_bank() { return nu_bank_; }
NuBank& nu_bank(int i) { return nu_bank_[i]; }
// Interim pulse height tally storage
vector<double>& pht_storage() { return pht_storage_; }
// Global tally accumulators
double& keff_tally_absorption() { return keff_tally_absorption_; }
double& keff_tally_collision() { return keff_tally_collision_; }
double& keff_tally_tracklength() { return keff_tally_tracklength_; }
double& keff_tally_leakage() { return keff_tally_leakage_; }
// Shows debug info
bool& trace() { return trace_; }
// Distance to the next collision
double& collision_distance() { return collision_distance_; }
// Number of events particle has undergone
int& n_event() { return n_event_; }
// Number of times variance reduction has caused a particle split
int n_split() const { return n_split_; }
int& n_split() { return n_split_; }
// Particle-specific factor for on-the-fly weight window adjustment
double ww_factor() const { return ww_factor_; }
double& ww_factor() { return ww_factor_; }
#ifdef DAGMC
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
#endif
// Number of progeny produced by this particle
int64_t& n_progeny() { return n_progeny_; }
// 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; }
//! Gets the pointer to the particle's current PRN seed
uint64_t* current_seed() { return seeds_ + stream_; }
const uint64_t* current_seed() const { return seeds_ + stream_; }
@ -513,14 +619,6 @@ public:
micro.last_E = 0.0;
}
//! resets all coordinate levels for the particle
void clear()
{
for (auto& level : coord_)
level.reset();
n_coord_ = 1;
}
//! Get track information based on particle's current state
TrackState get_track_state() const;

8
include/openmc/plot.h
View file

@ -121,7 +121,7 @@ struct IdData {
IdData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
// Members
@ -133,7 +133,7 @@ struct PropertyData {
PropertyData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
// Members
@ -205,7 +205,7 @@ T SlicePlotBase::get_map() const
#pragma omp parallel
{
Particle p;
GeometryState p;
p.r() = xyz;
p.u() = dir;
p.coord(0).universe = model::root_universe;
@ -291,7 +291,7 @@ private:
* find a distance to the boundary in a non-standard surface intersection
* check. It's an exhaustive search over surfaces in the top-level universe.
*/
static int advance_to_boundary_from_void(Particle& p);
static int advance_to_boundary_from_void(GeometryState& p);
/* Checks if a vector of two TrackSegments is equivalent. We define this
* to mean not having matching intersection lengths, but rather having

7
include/openmc/surface.h
View file

@ -105,12 +105,13 @@ public:
//! Determine the direction of a ray reflected from the surface.
//! \param[in] r The point at which the ray is incident.
//! \param[in] u Incident direction of the ray
//! \param[inout] p Pointer to the particle
//! \param[inout] p Pointer to the particle. Only DAGMC uses this.
//! \return Outgoing direction of the ray
virtual Direction reflect(Position r, Direction u, Particle* p) const;
virtual Direction reflect(
Position r, Direction u, GeometryState* p = nullptr) const;
virtual Direction diffuse_reflect(
Position r, Direction u, uint64_t* seed) const;
Position r, Direction u, uint64_t* seed, GeometryState* p = nullptr) const;
//! Evaluate the equation describing the surface.
//!

3
include/openmc/universe.h
View file

@ -9,6 +9,7 @@ namespace openmc {
class DAGUniverse;
#endif
class GeometryState;
class Universe;
class UniversePartitioner;
@ -32,7 +33,7 @@ public:
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
virtual bool find_cell(GeometryState& p) const;
BoundingBox bounding_box() const;

7
src/cell.cpp
View file

@ -784,7 +784,7 @@ std::string Region::str() const
//==============================================================================
std::pair<double, int32_t> Region::distance(
Position r, Direction u, int32_t on_surface, Particle* p) const
Position r, Direction u, int32_t on_surface) const
{
double min_dist {INFTY};
int32_t i_surf {std::numeric_limits<int32_t>::max()};
@ -1301,14 +1301,15 @@ vector<ParentCell> Cell::find_parent_cells(
{
// create a temporary particle
Particle dummy_particle {};
GeometryState dummy_particle {};
dummy_particle.r() = r;
dummy_particle.u() = {0., 0., 1.};
return find_parent_cells(instance, dummy_particle);
}
vector<ParentCell> Cell::find_parent_cells(int32_t instance, Particle& p) const
vector<ParentCell> Cell::find_parent_cells(
int32_t instance, GeometryState& p) const
{
// look up the particle's location
exhaustive_find_cell(p);

12
src/dagmc.cpp
View file

@ -404,7 +404,7 @@ int32_t DAGUniverse::implicit_complement_idx() const
return cell_idx_offset_ + dagmc_instance_->index_by_handle(ic) - 1;
}
bool DAGUniverse::find_cell(Particle& p) const
bool DAGUniverse::find_cell(GeometryState& p) const
{
// if the particle isn't in any of the other DagMC
// cells, place it in the implicit complement
@ -583,9 +583,8 @@ DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
};
std::pair<double, int32_t> DAGCell::distance(
Position r, Direction u, int32_t on_surface, Particle* p) const
Position r, Direction u, int32_t on_surface, GeometryState* p) const
{
Expects(p);
// if we've changed direction or we're not on a surface,
// reset the history and update last direction
if (u != p->last_dir()) {
@ -635,10 +634,9 @@ std::pair<double, int32_t> DAGCell::distance(
p->material() == MATERIAL_VOID
? "-1 (VOID)"
: std::to_string(model::materials[p->material()]->id());
auto lost_particle_msg = fmt::format(
p->mark_as_lost(fmt::format(
"No intersection found with DAGMC cell {}, filled with material {}", id_,
material_id);
p->mark_as_lost(lost_particle_msg);
material_id));
}
return {dist, surf_idx};
@ -723,7 +721,7 @@ Direction DAGSurface::normal(Position r) const
return dir;
}
Direction DAGSurface::reflect(Position r, Direction u, Particle* p) const
Direction DAGSurface::reflect(Position r, Direction u, GeometryState* p) const
{
Expects(p);
p->history().reset_to_last_intersection();

67
src/geometry.cpp
View file

@ -32,7 +32,7 @@ vector<int64_t> overlap_check_count;
// Non-member functions
//==============================================================================
bool check_cell_overlap(Particle& p, bool error)
bool check_cell_overlap(GeometryState& p, bool error)
{
int n_coord = p.n_coord();
@ -63,7 +63,7 @@ bool check_cell_overlap(Particle& p, bool error)
//==============================================================================
int cell_instance_at_level(const Particle& p, int level)
int cell_instance_at_level(const GeometryState& p, int level)
{
// throw error if the requested level is too deep for the geometry
if (level > model::n_coord_levels) {
@ -99,7 +99,8 @@ int cell_instance_at_level(const Particle& p, int level)
//==============================================================================
bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
bool find_cell_inner(
GeometryState& p, const NeighborList* neighbor_list, bool verbose)
{
// Find which cell of this universe the particle is in. Use the neighbor list
// to shorten the search if one was provided.
@ -156,7 +157,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
i_cell = p.lowest_coord().cell;
// Announce the cell that the particle is entering.
if (found && (settings::verbosity >= 10 || p.trace())) {
if (found && verbose) {
auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_);
write_message(msg, 1);
}
@ -243,10 +244,9 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
if (lat.outer_ != NO_OUTER_UNIVERSE) {
coord.universe = lat.outer_;
} else {
warning(fmt::format("Particle {} is outside lattice {} but the "
"lattice has no defined outer universe.",
p.mark_as_lost(fmt::format(
"Particle {} left lattice {}, but it has no outer definition.",
p.id(), lat.id_));
return false;
}
}
}
@ -259,7 +259,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
//==============================================================================
bool neighbor_list_find_cell(Particle& p)
bool neighbor_list_find_cell(GeometryState& p, bool verbose)
{
// Reset all the deeper coordinate levels.
@ -274,20 +274,20 @@ bool neighbor_list_find_cell(Particle& p)
// Search for the particle in that cell's neighbor list. Return if we
// found the particle.
bool found = find_cell_inner(p, &c.neighbors_);
bool found = find_cell_inner(p, &c.neighbors_, verbose);
if (found)
return found;
// The particle could not be found in the neighbor list. Try searching all
// cells in this universe, and update the neighbor list if we find a new
// neighboring cell.
found = find_cell_inner(p, nullptr);
found = find_cell_inner(p, nullptr, verbose);
if (found)
c.neighbors_.push_back(p.coord(coord_lvl).cell);
return found;
}
bool exhaustive_find_cell(Particle& p)
bool exhaustive_find_cell(GeometryState& p, bool verbose)
{
int i_universe = p.lowest_coord().universe;
if (i_universe == C_NONE) {
@ -299,17 +299,17 @@ bool exhaustive_find_cell(Particle& p)
for (int i = p.n_coord(); i < model::n_coord_levels; i++) {
p.coord(i).reset();
}
return find_cell_inner(p, nullptr);
return find_cell_inner(p, nullptr, verbose);
}
//==============================================================================
void cross_lattice(Particle& p, const BoundaryInfo& boundary)
void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
{
auto& coord {p.lowest_coord()};
auto& lat {*model::lattices[coord.lattice]};
if (settings::verbosity >= 10 || p.trace()) {
if (verbose) {
write_message(
fmt::format(" Crossing lattice {}. Current position ({},{},{}). r={}",
lat.id_, coord.lattice_i[0], coord.lattice_i[1], coord.lattice_i[2],
@ -336,10 +336,11 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary)
// The particle is outside the lattice. Search for it from the base coords.
p.n_coord() = 1;
bool found = exhaustive_find_cell(p);
if (!found && p.alive()) {
p.mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary",
p.id()));
if (!found) {
p.mark_as_lost(fmt::format("Particle {} could not be located after "
"crossing a boundary of lattice {}",
p.id(), lat.id_));
}
} else {
@ -352,10 +353,10 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary)
// this case, search for it from the base coords.
p.n_coord() = 1;
bool found = exhaustive_find_cell(p);
if (!found && p.alive()) {
p.mark_as_lost(fmt::format("Could not locate particle {} after "
"crossing a lattice boundary",
p.id()));
if (!found) {
p.mark_as_lost(fmt::format("Particle {} could not be located after "
"crossing a boundary of lattice {}",
p.id(), lat.id_));
}
}
}
@ -363,7 +364,7 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary)
//==============================================================================
BoundaryInfo distance_to_boundary(Particle& p)
BoundaryInfo distance_to_boundary(GeometryState& p)
{
BoundaryInfo info;
double d_lat = INFINITY;
@ -408,8 +409,9 @@ BoundaryInfo distance_to_boundary(Particle& p)
level_lat_trans = lattice_distance.second;
if (d_lat < 0) {
p.mark_as_lost(fmt::format(
"Particle {} had a negative distance to a lattice boundary", p.id()));
p.mark_as_lost(fmt::format("Particle {} had a negative distance "
"to a lattice boundary.",
p.id()));
}
}
@ -463,18 +465,19 @@ BoundaryInfo distance_to_boundary(Particle& p)
extern "C" int openmc_find_cell(
const double* xyz, int32_t* index, int32_t* instance)
{
Particle p;
GeometryState geom_state;
p.r() = Position {xyz};
p.u() = {0.0, 0.0, 1.0};
geom_state.r() = Position {xyz};
geom_state.u() = {0.0, 0.0, 1.0};
if (!exhaustive_find_cell(p)) {
set_errmsg(fmt::format("Could not find cell at position {}.", p.r()));
if (!exhaustive_find_cell(geom_state)) {
set_errmsg(
fmt::format("Could not find cell at position {}.", geom_state.r()));
return OPENMC_E_GEOMETRY;
}
*index = p.lowest_coord().cell;
*instance = p.cell_instance();
*index = geom_state.lowest_coord().cell;
*instance = geom_state.cell_instance();
return 0;
}

18
src/particle.cpp
View file

@ -277,7 +277,9 @@ void Particle::event_cross_surface()
boundary().lattice_translation[1] != 0 ||
boundary().lattice_translation[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(*this, boundary());
bool verbose = settings::verbosity >= 10 || trace();
cross_lattice(*this, boundary(), verbose);
event() = TallyEvent::LATTICE;
} else {
// Particle crosses surface
@ -406,7 +408,8 @@ void Particle::event_revive_from_secondary()
// have to determine it before the energy of the secondary particle can be
// removed from the pulse-height of this cell.
if (lowest_coord().cell == C_NONE) {
if (!exhaustive_find_cell(*this)) {
bool verbose = settings::verbosity >= 10 || trace();
if (!exhaustive_find_cell(*this, verbose)) {
mark_as_lost("Could not find the cell containing particle " +
std::to_string(id()));
return;
@ -556,7 +559,8 @@ void Particle::cross_surface()
}
#endif
if (neighbor_list_find_cell(*this))
bool verbose = settings::verbosity >= 10 || trace();
if (neighbor_list_find_cell(*this, verbose))
return;
// ==========================================================================
@ -564,7 +568,7 @@ void Particle::cross_surface()
// Remove lower coordinate levels
n_coord() = 1;
bool found = exhaustive_find_cell(*this);
bool found = exhaustive_find_cell(*this, verbose);
if (settings::run_mode != RunMode::PLOTTING && (!found)) {
// If a cell is still not found, there are two possible causes: 1) there is
@ -579,7 +583,7 @@ void Particle::cross_surface()
// Couldn't find next cell anywhere! This probably means there is an actual
// undefined region in the geometry.
if (!exhaustive_find_cell(*this)) {
if (!exhaustive_find_cell(*this, verbose)) {
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.");
@ -654,8 +658,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
// (unless we're using a dagmc model, which has exactly one universe)
n_coord() = 1;
if (surf.geom_type_ != GeometryType::DAG && !neighbor_list_find_cell(*this)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface " +
std::to_string(surf.id_) + ".");
mark_as_lost("Couldn't find particle after reflecting from surface " +
std::to_string(surf.id_) + ".");
return;
}

23
src/particle_data.cpp
View file

@ -1,6 +1,9 @@
#include "openmc/particle_data.h"
#include <sstream>
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
#include "openmc/nuclide.h"
@ -12,6 +15,21 @@
namespace openmc {
void GeometryState::mark_as_lost(const std::string& message)
{
mark_as_lost(message.c_str());
}
void GeometryState::mark_as_lost(const std::stringstream& message)
{
mark_as_lost(message.str());
}
void GeometryState::mark_as_lost(const char* message)
{
fatal_error(message);
}
void LocalCoord::rotate(const vector<double>& rotation)
{
r = r.rotate(rotation);
@ -30,13 +48,16 @@ void LocalCoord::reset()
rotated = false;
}
ParticleData::ParticleData()
GeometryState::GeometryState()
{
// Create and clear coordinate levels
coord_.resize(model::n_coord_levels);
cell_last_.resize(model::n_coord_levels);
clear();
}
ParticleData::ParticleData()
{
zero_delayed_bank();
// Every particle starts with no accumulated flux derivative. Note that in

43
src/plot.cpp
View file

@ -45,7 +45,7 @@ constexpr int32_t OVERLAP {-3};
IdData::IdData(size_t h_res, size_t v_res) : data_({v_res, h_res, 3}, NOT_FOUND)
{}
void IdData::set_value(size_t y, size_t x, const Particle& p, int level)
void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level)
{
// set cell data
if (p.n_coord() <= level) {
@ -78,7 +78,8 @@ PropertyData::PropertyData(size_t h_res, size_t v_res)
: data_({v_res, h_res, 2}, NOT_FOUND)
{}
void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level)
void PropertyData::set_value(
size_t y, size_t x, const GeometryState& p, int level)
{
Cell* c = model::cells.at(p.lowest_coord().cell).get();
data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
@ -1084,7 +1085,7 @@ void ProjectionPlot::set_output_path(pugi::xml_node node)
// Advances to the next boundary from outside the geometry
// Returns -1 if no intersection found, and the surface index
// if an intersection was found.
int ProjectionPlot::advance_to_boundary_from_void(Particle& p)
int ProjectionPlot::advance_to_boundary_from_void(GeometryState& p)
{
constexpr double scoot = 1e-5;
double min_dist = {INFINITY};
@ -1234,20 +1235,8 @@ void ProjectionPlot::create_output() const
const int n_threads = num_threads();
const int tid = thread_num();
SourceSite s; // Where particle starts from (camera)
s.E = 1;
s.wgt = 1;
s.delayed_group = 0;
s.particle = ParticleType::photon; // just has to be something reasonable
s.parent_id = 1;
s.progeny_id = 2;
s.r = camera_position_;
Particle p;
s.u.x = 1.0;
s.u.y = 0.0;
s.u.z = 0.0;
p.from_source(&s);
GeometryState p;
p.u() = {1.0, 0.0, 0.0};
int vert = tid;
for (int iter = 0; iter <= pixels_[1] / n_threads; iter++) {
@ -1263,6 +1252,10 @@ void ProjectionPlot::create_output() const
for (int horiz = 0; horiz < pixels_[0]; ++horiz) {
// Projection mode below decides ray starting conditions
Position init_r;
Direction init_u;
// Generate the starting position/direction of the ray
if (orthographic_width_ == 0.0) { // perspective projection
double this_phi =
@ -1273,18 +1266,22 @@ void ProjectionPlot::create_output() const
camera_local_vec.x = std::cos(this_phi) * std::sin(this_mu);
camera_local_vec.y = std::sin(this_phi) * std::sin(this_mu);
camera_local_vec.z = std::cos(this_mu);
s.u = camera_local_vec.rotate(camera_to_model);
init_u = camera_local_vec.rotate(camera_to_model);
init_r = camera_position_;
} else { // orthographic projection
s.u = looking_direction;
init_u = looking_direction;
double x_pix_coord = (static_cast<double>(horiz) - p0 / 2.0) / p0;
double y_pix_coord = (static_cast<double>(vert) - p1 / 2.0) / p0;
s.r = camera_position_ +
cam_yaxis * x_pix_coord * orthographic_width_ +
cam_zaxis * y_pix_coord * orthographic_width_;
init_r = camera_position_;
init_r += cam_yaxis * x_pix_coord * orthographic_width_;
init_r += cam_zaxis * y_pix_coord * orthographic_width_;
}
p.from_source(&s); // put particle at camera
// Resets internal geometry state of particle
p.init_from_r_u(init_r, init_u);
bool hitsomething = false;
bool intersection_found = true;
int loop_counter = 0;

4
src/surface.cpp
View file

@ -129,7 +129,7 @@ bool Surface::sense(Position r, Direction u) const
return f > 0.0;
}
Direction Surface::reflect(Position r, Direction u, Particle* p) const
Direction Surface::reflect(Position r, Direction u, GeometryState* p) const
{
// Determine projection of direction onto normal and squared magnitude of
// normal.
@ -140,7 +140,7 @@ Direction Surface::reflect(Position r, Direction u, Particle* p) const
}
Direction Surface::diffuse_reflect(
Position r, Direction u, uint64_t* seed) const
Position r, Direction u, uint64_t* seed, GeometryState* p) const
{
// Diffuse reflect direction according to the normal.
// cosine distribution

3
src/universe.cpp
View file

@ -3,6 +3,7 @@
#include <set>
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
namespace openmc {
@ -36,7 +37,7 @@ void Universe::to_hdf5(hid_t universes_group) const
close_group(group);
}
bool Universe::find_cell(Particle& p) const
bool Universe::find_cell(GeometryState& p) const
{
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};