openmc-dev / openmc / 18542382223

#ifndef OPENMC_PARTICLE_DATA_H

#define OPENMC_PARTICLE_DATA_H

#include "openmc/array.h"

#include "openmc/constants.h"

#include "openmc/position.h"

#include "openmc/random_lcg.h"

#include "openmc/tallies/filter_match.h"

#include "openmc/vector.h"

#ifdef OPENMC_DAGMC_ENABLED

#include "DagMC.hpp"

#endif

namespace openmc {

//==============================================================================

// Constants

//==============================================================================

// Since cross section libraries come with different numbers of delayed groups

// (e.g. ENDF/B-VII.1 has 6 and JEFF 3.1.1 has 8 delayed groups) and we don't

// yet know what cross section library is being used when the tallies.xml file

// is read in, we want to have an upper bound on the size of the array we

// use to store the bins for delayed group tallies.

constexpr int MAX_DELAYED_GROUPS {8};

constexpr double CACHE_INVALID {-1.0};

//==========================================================================

// Aliases and type definitions

//! Particle types

enum class ParticleType { neutron, photon, electron, positron };

//! Saved ("banked") state of a particle

//! NOTE: This structure's MPI type is built in initialize_mpi() of

//! initialize.cpp. Any changes made to the struct here must also be

//! made when building the Bank MPI type in initialize_mpi().

//! NOTE: This structure is also used on the python side, and is defined

//! in lib/core.py. Changes made to the type here must also be made to the

//! python defintion.

struct SourceSite {

  Position r;

  Direction u;

  double E;

  double time {0.0};

  double wgt {1.0};

  int delayed_group {0};

  int surf_id {SURFACE_NONE};

  ParticleType particle;

  // Extra attributes that don't show up in source written to file

  int parent_nuclide {-1};

  int64_t parent_id;

  int64_t progeny_id;

};

//! State of a particle used for particle track files

struct TrackState {

  Position r;           //!< Position in [cm]

  Direction u;          //!< Direction

  double E;             //!< Energy in [eV]

  double time {0.0};    //!< Time in [s]

  double wgt {1.0};     //!< Weight

  int cell_id;          //!< Cell ID

  int cell_instance;    //!< Cell instance

  int material_id {-1}; //!< Material ID (default value indicates void)

};

//! Full history of a single particle's track states

struct TrackStateHistory {

  ParticleType particle;

  std::vector<TrackState> states;

};

//! Saved ("banked") state of a particle, for nu-fission tallying

struct NuBank {

  double E;          //!< particle energy

  double wgt;        //!< particle weight

  int delayed_group; //!< particle delayed group

};

class LocalCoord {

public:

  void rotate(const vector<double>& rotation);

  //! clear data from a single coordinate level

  void reset();

  // accessors

private:

  // Data members

  Position r_;  //!< particle position

  Direction u_; //!< particle direction

  int cell_ {-1};

  int universe_ {-1};

  int lattice_ {-1};

  array<int, 3> lattice_index_ {{-1, -1, -1}};

  bool rotated_ {false}; //!< Is the level rotated?

};

//==============================================================================

//! Cached microscopic cross sections for a particular nuclide at the current

//! energy

//==============================================================================

struct NuclideMicroXS {

  // Microscopic cross sections in barns

  double total;      //!< total cross section

  double absorption; //!< absorption (disappearance)

  double fission;    //!< fission

  double nu_fission; //!< neutron production from fission

  double elastic;         //!< If sab_frac is not 1 or 0, then this value is

                          //!<   averaged over bound and non-bound nuclei

  double thermal;         //!< Bound thermal elastic & inelastic scattering

  double thermal_elastic; //!< Bound thermal elastic scattering

  double photon_prod;     //!< microscopic photon production xs

  // Cross sections for depletion reactions (note that these are not stored in

  // macroscopic cache)

  double reaction[DEPLETION_RX.size()];

  // Indicies and factors needed to compute cross sections from the data tables

  int index_grid;       //!< Index on nuclide energy grid

  int index_temp;       //!< Temperature index for nuclide

  double interp_factor; //!< Interpolation factor on nuc. energy grid

  int index_sab {-1};   //!< Index in sab_tables

  int index_temp_sab;   //!< Temperature index for sab_tables

  double sab_frac;      //!< Fraction of atoms affected by S(a,b)

  bool use_ptable;      //!< In URR range with probability tables?

  // Energy and temperature last used to evaluate these cross sections.  If

  // these values have changed, then the cross sections must be re-evaluated.

  double last_E {0.0};       //!< Last evaluated energy

  double last_sqrtkT {0.0};  //!< Last temperature in sqrt(Boltzmann constant

                             //!< * temperature (eV))

  double ncrystal_xs {-1.0}; //!< NCrystal cross section

};

//==============================================================================

//! Cached microscopic photon cross sections for a particular element at the

//! current energy

//==============================================================================

struct ElementMicroXS {

  int index_grid;         //!< index on element energy grid

  double last_E {0.0};    //!< last evaluated energy in [eV]

  double interp_factor;   //!< interpolation factor on energy grid

  double total;           //!< microscopic total photon xs

  double coherent;        //!< microscopic coherent xs

  double incoherent;      //!< microscopic incoherent xs

  double photoelectric;   //!< microscopic photoelectric xs

  double pair_production; //!< microscopic pair production xs

};

//==============================================================================

// MacroXS contains cached macroscopic cross sections for the material a

// particle is traveling through

//==============================================================================

struct MacroXS {

  double total;       //!< macroscopic total xs

  double absorption;  //!< macroscopic absorption xs

  double fission;     //!< macroscopic fission xs

  double nu_fission;  //!< macroscopic production xs

  double photon_prod; //!< macroscopic photon production xs

  // Photon cross sections

  double coherent;        //!< macroscopic coherent xs

  double incoherent;      //!< macroscopic incoherent xs

  double photoelectric;   //!< macroscopic photoelectric xs

  double pair_production; //!< macroscopic pair production xs

};

//==============================================================================

// Cache contains the cached data for an MGXS object

//==============================================================================

struct CacheDataMG {

  int material {-1}; //!< material index

  double sqrtkT;     //!< last temperature corresponding to t

  int t {0};         //!< temperature index

  int a {0};         //!< angle index

  Direction u;       //!< angle that corresponds to a

};

//==============================================================================

// Information about nearest boundary crossing

//==============================================================================

class BoundaryInfo {

public:

  {

  const double& distance() const { return distance_; }

  int coord_level() const { return coord_level_; }

  {

  }

  // TODO: off-by-one

private:

  // Data members

  double distance_ {INFINITY}; //!< distance to nearest boundary

  int surface_ {

    SURFACE_NONE};      //!< surface token, non-zero if boundary is surface

  int coord_level_ {0}; //!< coordinate level after crossing boundary

  array<int, 3> lattice_translation_ {

    0, 0, 0}; //!< 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

  {

    }

    }

  //! moves the particle by the specified distance to its next location

  //! \param distance the distance the particle is moved

  void move_distance(double distance);

  void advance_to_boundary_from_void();

  // Initialize all internal state from position and direction

  {

  // 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.

  // Number of current coordinate levels

  // Offset for distributed properties

  // Coordinates for all nesting levels

  // Innermost universe nesting coordinates

  // Last coordinates on all nesting levels, before crossing a surface

  // Coordinates at birth

  // Coordinates of last collision or reflective/periodic surface

  // crossing for current tallies

  // Previous direction and spatial coordinates before a collision

  // Accessors for position in global coordinates

  // Accessors for position in local coordinates

  const Position& r_local() const { return coord_[n_coord_ - 1].r(); }

  // Accessors for direction in global coordinates

  // Accessors for direction in local coordinates

  const Direction& u_local() const { return coord_[n_coord_ - 1].u(); }

  // Surface token for the surface that the particle is currently on

  // Surface index based on the current value of the surface_ attribute

  {

    // TODO: off-by-one

  }

  // Boundary information

#ifdef OPENMC_DAGMC_ENABLED

  // DagMC state variables

#endif

  // material of current and last cell

  // temperature of current and last cell

  // density multiplier of the current and last cell

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_born_;         //!< coordinates at birth

  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_ {

    SURFACE_NONE}; //!< surface token for surface the particle is currently 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

  double density_mult_ {1.0};      //!< density multiplier

  double density_mult_last_ {1.0}; //!< last density multiplier

#ifdef OPENMC_DAGMC_ENABLED

  moab::DagMC::RayHistory history_;

  Direction last_dir_;

#endif

};

//============================================================================

//! Defines how particle data is laid out in memory

//============================================================================

/*

 * This class was added in order to separate the layout and access of particle

 * data from particle physics operations during a development effort to get

 * OpenMC running on GPUs. In the event-based Monte Carlo method, one creates

 * an array of particles on which actions like cross section lookup and surface

 * crossing are done en masse, which works best on vector computers of yore and

 * modern GPUs. It has been shown in the below publication [1] that arranging

 * particle data into a structure of arrays rather than an array of structures

 * enhances performance on GPUs. For instance, rather than having an

 * std::vector<Particle> where consecutive particle energies would be separated

 * by about 400 bytes, one would create a structure which has a single

 * std::vector<double> of energies.  The motivation here is that more coalesced

 * memory accesses occur, in the parlance of GPU programming.

 *

 * So, this class enables switching between the array-of-structures and

 * structure- of-array data layout at compile time. In GPU branches of the

 * code, our Particle class inherits from a class that provides an array of

 * particle energies, and can access them using the E() method (defined below).

 * In the CPU code, we inherit from this class which gives the conventional

 * layout of particle data, useful for history-based tracking.

 *

 * As a result, we always use the E(), r_last(), etc. methods to access

 * particle data in order to keep a unified interface between

 * structure-of-array and array-of-structure code on either CPU or GPU code

 * while sharing the same physics code on each codebase.

 *

 * [1] Hamilton, Steven P., Stuart R. Slattery, and Thomas M. Evans.

 *   “Multigroup Monte Carlo on GPUs: Comparison of History- and Event-Based

 *   Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.

 *   https://doi.org/10.1016/j.anucene.2017.11.032.

 */

class ParticleData : public GeometryState {

private:

  //==========================================================================

  // Data members -- see public: below for descriptions

  vector<NuclideMicroXS> neutron_xs_;

  vector<ElementMicroXS> photon_xs_;

  MacroXS macro_xs_;

  CacheDataMG mg_xs_cache_;

  ParticleType type_ {ParticleType::neutron};

  double E_;

  double E_last_;

  int g_ {0};

  int g_last_;

  double wgt_ {1.0};

  double wgt_born_ {1.0};

  double wgt_ww_born_ {-1.0};

  double mu_;

  double time_ {0.0};

  double time_last_ {0.0};

  double wgt_last_ {1.0};

  bool fission_ {false};

  TallyEvent event_;

  int event_nuclide_;

  int event_mt_;

  int delayed_group_ {0};

  int parent_nuclide_ {-1};

  int n_bank_ {0};

  double bank_second_E_ {0.0};

  double wgt_bank_ {0.0};

  int n_delayed_bank_[MAX_DELAYED_GROUPS];

  int cell_born_ {-1};

  // Iterated Fission Probability

  double lifetime_ {0.0}; //!< neutron lifetime [s]

  int n_collision_ {0};

  bool write_track_ {false};

  uint64_t seeds_[N_STREAMS];

  int stream_;

  vector<SourceSite> secondary_bank_;

  int64_t current_work_;

  vector<double> flux_derivs_;

  vector<FilterMatch> filter_matches_;

  vector<TrackStateHistory> tracks_;

  vector<NuBank> nu_bank_;

  vector<double> pht_storage_;

  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};

  double collision_distance_;

  int n_event_ {0};

  int n_split_ {0};

  double ww_factor_ {0.0};

  int64_t n_progeny_ {0};

public:

  //----------------------------------------------------------------------------

  // Constructors

  ParticleData();

  //==========================================================================

  // Methods and accessors

  // Cross section caches

  {

  } // Microscopic neutron cross sections

  // Microscopic photon cross sections

  // Macroscopic cross sections

  // Multigroup macroscopic cross sections

  // Particle type (n, p, e, gamma, etc)

  // Current particle energy, energy before collision,

  // and corresponding multigroup group indices. Energy

  // units are eV.

  // Statistic weight of particle. Setting to zero indicates that the particle

  // is dead.

  // Statistic weight of particle at birth

  double wgt_born() const { return wgt_born_; }

  // Weight window value at birth

  const double& wgt_ww_born() const { return wgt_ww_born_; }

  // Statistic weight of particle at last collision

  // Whether particle is alive

  // Polar scattering angle after a collision

  // Tracks the time of a particle as it traverses the problem.

  // Units are seconds.

  // Particle lifetime

  // What event took place, described in greater detail below

  // Post-collision data

  {

  } // energy of last reaction secondaries

  {

  } // number of delayed fission sites

  {

  } // number of delayed fission sites

  // Index of cell particle is born in

  // Total number of collisions suffered by particle

  // whether this track is to be written

  // RNG state

  // secondary particle bank

  SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }

  // Current simulation work index

  // Used in tally derivatives

  // Matches of tallies

  // Tracks to output to file

  // Bank of recently fissioned particles

  // Interim pulse height tally storage

  // Global tally accumulators

  // Shows debug info

  // Distance to the next collision

  // Number of events particle has undergone

  // Number of times variance reduction has caused a particle split

  int n_split() const { return n_split_; }

  // Particle-specific factor for on-the-fly weight window adjustment

  double ww_factor() const { return ww_factor_; }

  // Number of progeny produced by this particle

  //! Gets the pointer to the particle's current PRN seed

  const uint64_t* current_seed() const { return seeds_ + stream_; }

  //! Force recalculation of neutron xs by setting last energy to zero

  {

  //! Get track information based on particle's current state

  TrackState get_track_state() const;

  {

    }

  {

    }

};

} // namespace openmc

#endif // OPENMC_PARTICLE_DATA_H