crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1183 from paulromano/cpp-cleanup

Browse source 
General improvements in C++ code
This commit is contained in:
Sterling Harper 2019-03-05 09:39:54 -05:00 • committed by GitHub
commit e6ae490a31
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
78 changed files with 1356 additions and 1422 deletions
Download patch file Download diff file Expand all files Collapse all files

22
include/openmc/bank.h
View file

@ -4,23 +4,13 @@
#include <cstdint>
#include <vector>
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<Bank>, 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<openmc::Bank>;
extern template class std::vector<openmc::Particle::Bank>;
namespace openmc {
@ -32,10 +22,10 @@ namespace simulation {
extern "C" int64_t n_bank;
extern std::vector<Bank> source_bank;
extern std::vector<Bank> fission_bank;
extern std::vector<Particle::Bank> source_bank;
extern std::vector<Particle::Bank> fission_bank;
#ifdef _OPENMP
extern std::vector<Bank> master_fission_bank;
extern std::vector<Particle::Bank> master_fission_bank;
#endif
#pragma omp threadprivate(fission_bank, n_bank)

18
include/openmc/capi.h
View file

@ -6,23 +6,7 @@
#include <stddef.h>
#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);

13
include/openmc/cell.h
View file

@ -3,6 +3,7 @@
#include <cstdint>
#include <limits>
#include <memory> // for unique_ptr
#include <string>
#include <unordered_map>
#include <vector>
@ -44,15 +45,11 @@ class Cell;
class Universe;
namespace model {
extern std::vector<std::unique_ptr<Cell>> cells;
extern std::unordered_map<int32_t, int32_t> cell_map;
extern "C" int32_t n_cells;
extern std::vector<Cell*> cells;
extern std::unordered_map<int32_t, int32_t> cell_map;
extern std::vector<Universe*> universes;
extern std::unordered_map<int32_t, int32_t> universe_map;
extern std::vector<std::unique_ptr<Universe>> universes;
extern std::unordered_map<int32_t, int32_t> universe_map;
} // namespace model
//==============================================================================

10
include/openmc/hdf5_interface.h
View file

@ -332,7 +332,17 @@ void read_dataset(hid_t obj_id, const char* name, xt::xtensor<T, N>& 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<double, 3> x;
read_dataset(obj_id, name, x, indep);
r.x = x[0];
r.y = x[1];
r.z = x[2];
}
template <typename T, std::size_t N>

13
include/openmc/lattice.h
View file

@ -3,6 +3,7 @@
#include <array>
#include <cstdint>
#include <memory> // for unique_ptr
#include <string>
#include <unordered_map>
#include <vector>
@ -33,10 +34,8 @@ enum class LatticeType {
class Lattice;
namespace model {
extern std::vector<Lattice*> lattices;
extern std::unordered_map<int32_t, int32_t> lattice_map;
extern std::vector<std::unique_ptr<Lattice>> lattices;
extern std::unordered_map<int32_t, int32_t> lattice_map;
} // namespace model
//==============================================================================
@ -105,7 +104,7 @@ public:
//! \brief Find the lattice tile indices for a given point.
//! \param r A 3D Cartesian coordinate.
//! \return An array containing the indices of a lattice tile.
virtual std::array<int, 3> get_indices(Position r) const = 0;
virtual std::array<int, 3> get_indices(Position r, Direction u) const = 0;
//! \brief Get coordinates local to a lattice tile.
//! \param r A 3D Cartesian coordinate.
@ -213,7 +212,7 @@ public:
std::pair<double, std::array<int, 3>>
distance(Position r, Direction u, const std::array<int, 3>& i_xyz) const;
std::array<int, 3> get_indices(Position r) const;
std::array<int, 3> get_indices(Position r, Direction u) const;
Position
get_local_position(Position r, const std::array<int, 3> i_xyz) const;
@ -253,7 +252,7 @@ public:
std::pair<double, std::array<int, 3>>
distance(Position r, Direction u, const std::array<int, 3>& i_xyz) const;
std::array<int, 3> get_indices(Position r) const;
std::array<int, 3> get_indices(Position r, Direction u) const;
Position
get_local_position(Position r, const std::array<int, 3> i_xyz) const;

2
include/openmc/material.h
View file

@ -23,7 +23,7 @@ class Material;
namespace model {
extern std::vector<Material*> materials;
extern std::vector<std::unique_ptr<Material>> materials;
extern std::unordered_map<int32_t, int32_t> material_map;
} // namespace model

4
include/openmc/math_functions.h
View file

@ -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

2
include/openmc/mesh.h
View file

@ -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<double> count_sites(int64_t n, const Bank* bank,
xt::xarray<double> count_sites(int64_t n, const Particle::Bank* bank,
int n_energy, const double* energies, bool* outside) const;
int id_ {-1}; //!< User-specified ID

8
include/openmc/mgxs.h
View file

@ -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

2
include/openmc/mgxs_interface.h
View file

@ -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

136
include/openmc/particle.h
View file

@ -9,7 +9,7 @@
#include <sstream>
#include <string>
#include "openmc/capi.h"
#include "openmc/position.h"
namespace openmc {
@ -33,20 +33,19 @@ constexpr int MAX_LOST_PARTICLES {10};
// Maximum number of lost particles, relative to the total number of particles
constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
//! Particle types
enum class ParticleType {
neutron, photon, electron, positron
};
//==============================================================================
// Class declarations
//==============================================================================
struct LocalCoord {
Position r; //!< particle position
Direction u; //!< particle direction
int cell {-1};
int universe {-1};
int lattice {-1};
int lattice_x {-1};
int lattice_y {-1};
int lattice_z {-1};
double xyz[3]; //!< particle position
double uvw[3]; //!< particle direction
bool rotated {false}; //!< Is the level rotated?
//! clear data from a single coordinate level
@ -57,70 +56,105 @@ struct LocalCoord {
//! State of a particle being transported through geometry
//============================================================================
struct Particle {
int64_t id; //!< Unique ID
int type; //!< Particle type (n, p, e, etc.)
class Particle {
public:
//! Particle types
enum class Type {
neutron, photon, electron, positron
};
int n_coord; //!< number of current coordinate levels
int cell_instance; //!< offset for distributed properties
LocalCoord coord[MAX_COORD]; //!< coordinates for all levels
//! Saved ("banked") state of a particle
struct Bank {
Position r;
Direction u;
double E;
double wgt;
int delayed_group;
Type particle;
};
// Constructors
Particle();
int64_t id_; //!< Unique ID
Type type_ {Type::neutron}; //!< Particle type (n, p, e, etc.)
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels
// Particle coordinates before crossing a surface
int last_n_coord; //!< number of current coordinates
int last_cell[MAX_COORD]; //!< coordinates for all levels
int n_coord_last_ {1}; //!< number of current coordinates
int cell_last_[MAX_COORD]; //!< coordinates for all levels
// Energy data
double E; //!< post-collision energy in eV
double last_E; //!< pre-collision energy in eV
int g; //!< post-collision energy group (MG only)
int last_g; //!< pre-collision energy group (MG only)
double E_; //!< post-collision energy in eV
double E_last_; //!< pre-collision energy in eV
int g_ {0}; //!< post-collision energy group (MG only)
int g_last_; //!< pre-collision energy group (MG only)
// Other physical data
double wgt; //!< particle weight
double mu; //!< angle of scatter
bool alive; //!< is particle alive?
double wgt_ {1.0}; //!< particle weight
double mu_; //!< angle of scatter
bool alive_ {true}; //!< is particle alive?
// Other physical data
double last_xyz_current[3]; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
double last_xyz[3]; //!< previous coordinates
double last_uvw[3]; //!< previous direction coordinates
double last_wgt; //!< pre-collision particle weight
double absorb_wgt; //!< weight absorbed for survival biasing
Position r_last_current_; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
Position r_last_; //!< previous coordinates
Direction u_last_; //!< previous direction coordinates
double wgt_last_ {1.0}; //!< pre-collision particle weight
double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing
// What event took place
bool fission; //!< did particle cause implicit fission
int event; //!< scatter, absorption
int event_nuclide; //!< index in nuclides array
int event_MT; //!< reaction MT
int delayed_group; //!< delayed group
bool fission_ {false}; //!< did particle cause implicit fission
int event_; //!< scatter, absorption
int event_nuclide_; //!< index in nuclides array
int event_mt_; //!< reaction MT
int delayed_group_ {0}; //!< delayed group
// Post-collision physical data
int n_bank; //!< number of fission sites banked
double wgt_bank; //!< weight of fission sites banked
int n_delayed_bank[MAX_DELAYED_GROUPS]; //!< number of delayed fission
int n_bank_ {0}; //!< number of fission sites banked
double wgt_bank_ {0.0}; //!< weight of fission sites banked
int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission
//!< sites banked
// Indices for various arrays
int surface; //!< index for surface particle is on
int cell_born; //!< index for cell particle was born in
int material; //!< index for current material
int last_material; //!< index for last material
int surface_ {0}; //!< index for surface particle is on
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
// Temperature of current cell
double sqrtkT; //!< sqrt(k_Boltzmann * temperature) in eV
double last_sqrtkT; //!< last temperature
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
// Statistical data
int n_collision; //!< number of collisions
int n_collision_ {0}; //!< number of collisions
// Track output
bool write_track {false};
bool write_track_ {false};
// Secondary particles created
int64_t n_secondary {};
Bank secondary_bank[MAX_SECONDARY];
int64_t n_secondary_ {};
Bank secondary_bank_[MAX_SECONDARY];
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
//! resets all coordinate levels for the particle
void clear();
@ -129,14 +163,10 @@ struct Particle {
//
//! stores the current phase space attributes of the particle in the
//! secondary bank and increments the number of sites in the secondary bank.
//! \param uvw Direction of the secondary particle
//! \param u Direction of the secondary particle
//! \param E Energy of the secondary particle in [eV]
//! \param type Particle type
//! \param run_CE Whether continuous-energy data is being used
void create_secondary(const double* uvw, double E, int type, bool run_CE);
//! sets default attributes for a particle
void initialize();
void create_secondary(Direction u, double E, Type type);
//! initialize from a source site
//

14
include/openmc/physics.h
View file

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

2
include/openmc/physics_mg.h
View file

@ -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

5
include/openmc/position.h
View file

@ -1,6 +1,7 @@
#ifndef OPENMC_POSITION_H
#define OPENMC_POSITION_H
#include <array>
#include <cmath> // for sqrt
#include <stdexcept> // for out_of_range
#include <vector>
@ -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<double> xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
Position(const std::vector<double>& xyz) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { };
Position(const std::array<double, 3>& 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) {

2
include/openmc/reaction_product.h
View file

@ -44,7 +44,7 @@ public:
//! \param[out] mu Outgoing cosine with respect to current direction
void sample(double E_in, double& E_out, double& mu) const;
ParticleType particle_; //!< Particle type
Particle::Type particle_; //!< Particle type
EmissionMode emission_mode_; //!< Emission mode
double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]
std::unique_ptr<Function1D> yield_; //!< Yield as a function of energy

7
include/openmc/source.h
View file

@ -9,7 +9,6 @@
#include "pugixml.hpp"
#include "openmc/bank.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/particle.h"
@ -40,12 +39,12 @@ public:
//! Sample from the external source distribution
//! \return Sampled site
Bank sample() const;
Particle::Bank sample() const;
// Properties
double strength() const { return strength_; }
private:
ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
Particle::Type particle_ {Particle::Type::neutron}; //!< Type of particle emitted
double strength_ {1.0}; //!< Source strength
UPtrSpace space_; //!< Spatial distribution
UPtrAngle angle_; //!< Angular distribution
@ -62,7 +61,7 @@ extern "C" void initialize_source();
//! Sample a site from all external source distributions in proportion to their
//! source strength
//! \return Sampled source site
Bank sample_external_source();
Particle::Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();

9
include/openmc/surface.h
View file

@ -1,9 +1,10 @@
#ifndef OPENMC_SURFACE_H
#define OPENMC_SURFACE_H
#include <map>
#include <memory> // for unique_ptr
#include <limits> // For numeric_limits
#include <string>
#include <unordered_map>
#include <vector>
#include "hdf5.h"
@ -35,10 +36,8 @@ extern "C" const int BC_PERIODIC;
class Surface;
namespace model {
extern std::vector<Surface*> surfaces;
extern std::map<int, int> surface_map;
extern std::vector<std::unique_ptr<Surface>> surfaces;
extern std::unordered_map<int, int> surface_map;
} // namespace model
//==============================================================================

3
include/openmc/tallies/filter_particle.h
View file

@ -3,6 +3,7 @@
#include <vector>
#include "openmc/particle.h"
#include "openmc/tallies/filter.h"
namespace openmc {
@ -27,7 +28,7 @@ public:
std::string text_label(int bin) const override;
std::vector<int> particles_;
std::vector<Particle::Type> particles_;
};
} // namespace openmc

6
openmc/capi/core.py
View file

@ -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)]

8
openmc/cmfd.py
View file

@ -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

12
src/bank.cpp
View file

@ -6,7 +6,7 @@
#include <cstdint>
// Explicit template instantiation definition
template class std::vector<openmc::Bank>;
template class std::vector<openmc::Particle::Bank>;
namespace openmc {
@ -18,10 +18,10 @@ namespace simulation {
int64_t n_bank;
std::vector<Bank> source_bank;
std::vector<Bank> fission_bank;
std::vector<Particle::Bank> source_bank;
std::vector<Particle::Bank> fission_bank;
#ifdef _OPENMP
std::vector<Bank> master_fission_bank;
std::vector<Particle::Bank> 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.");

17
src/bremsstrahlung.cpp
View file

@ -28,20 +28,20 @@ std::vector<Bremsstrahlung> ttb;
void thick_target_bremsstrahlung(Particle& p, double* E_lost)
{
if (p.material == MATERIAL_VOID) return;
if (p.material_ == MATERIAL_VOID) return;
int photon = static_cast<int>(ParticleType::photon);
if (p.E < settings::energy_cutoff[photon]) return;
int photon = static_cast<int>(Particle::Type::photon);
if (p.E_ < settings::energy_cutoff[photon]) return;
// Get bremsstrahlung data for this material and particle type
BremsstrahlungData* mat;
if (p.type == static_cast<int>(ParticleType::positron)) {
mat = &model::materials[p.material]->ttb_->positron;
if (p.type_ == Particle::Type::positron) {
mat = &model::materials[p.material_]->ttb_->positron;
} else {
mat = &model::materials[p.material]->ttb_->electron;
mat = &model::materials[p.material_]->ttb_->electron;
}
double e = std::log(p.E);
double e = std::log(p.E_);
auto n_e = data::ttb_e_grid.size();
// Find the lower bounding index of the incident electron energy
@ -108,8 +108,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
if (w > settings::energy_cutoff[photon]) {
// Create secondary photon
int photon_ = static_cast<int>(ParticleType::photon);
p.create_secondary(p.coord[0].uvw, w, photon_, true);
p.create_secondary(p.u(), w, Particle::Type::photon);
*E_lost += w;
}
}

16
src/cell.cpp
View file

@ -22,13 +22,11 @@ namespace openmc {
//==============================================================================
namespace model {
std::vector<std::unique_ptr<Cell>> cells;
std::unordered_map<int32_t, int32_t> cell_map;
std::vector<Cell*> cells;
std::unordered_map<int32_t, int32_t> cell_map;
std::vector<Universe*> universes;
std::unordered_map<int32_t, int32_t> universe_map;
std::vector<std::unique_ptr<Universe>> universes;
std::unordered_map<int32_t, int32_t> universe_map;
} // namespace model
//==============================================================================
@ -653,7 +651,7 @@ void read_cells(pugi::xml_node node)
// Loop over XML cell elements and populate the array.
model::cells.reserve(n_cells);
for (pugi::xml_node cell_node : node.children("cell")) {
model::cells.push_back(new CSGCell(cell_node));
model::cells.push_back(std::make_unique<CSGCell>(cell_node));
}
// Fill the cell map.
@ -674,7 +672,7 @@ void read_cells(pugi::xml_node node)
int32_t uid = model::cells[i]->universe_;
auto it = model::universe_map.find(uid);
if (it == model::universe_map.end()) {
model::universes.push_back(new Universe());
model::universes.push_back(std::make_unique<Universe>());
model::universes.back()->id_ = uid;
model::universes.back()->cells_.push_back(i);
model::universe_map[uid] = model::universes.size() - 1;
@ -823,7 +821,7 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
if (index_start) *index_start = model::cells.size();
if (index_end) *index_end = model::cells.size() + n - 1;
for (int32_t i = 0; i < n; i++) {
model::cells.push_back(new CSGCell());
model::cells.push_back(std::make_unique<CSGCell>());
}
return 0;
}

10
src/cross_sections.cpp
View file

@ -228,7 +228,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
// Determine if minimum/maximum energy for this nuclide is greater/less
// than the previous
if (data::nuclides[i_nuclide]->grid_.size() >= 1) {
int neutron = static_cast<int>(ParticleType::neutron);
int neutron = static_cast<int>(Particle::Type::neutron);
data::energy_min[neutron] = std::max(data::energy_min[neutron],
data::nuclides[i_nuclide]->grid_[0].energy.front());
data::energy_max[neutron] = std::min(data::energy_max[neutron],
@ -262,7 +262,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
// the previous
const auto& elem {data::elements.back()};
if (elem.energy_.size() >= 1) {
int photon = static_cast<int>(ParticleType::photon);
int photon = static_cast<int>(Particle::Type::photon);
int n = elem.energy_.size();
data::energy_min[photon] = std::max(data::energy_min[photon],
std::exp(elem.energy_(1)));
@ -321,7 +321,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
for (auto& nuc : data::nuclides) {
nuc->init_grid();
}
int neutron = static_cast<int>(ParticleType::neutron);
int neutron = static_cast<int>(Particle::Type::neutron);
simulation::log_spacing = std::log(data::energy_max[neutron] /
data::energy_min[neutron]) / settings::n_log_bins;
@ -329,7 +329,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
// Determine if minimum/maximum energy for bremsstrahlung is greater/less
// than the current minimum/maximum
if (data::ttb_e_grid.size() >= 1) {
int photon = static_cast<int>(ParticleType::photon);
int photon = static_cast<int>(Particle::Type::photon);
int n_e = data::ttb_e_grid.size();
data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1));
data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1));
@ -345,7 +345,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
// grid has not been allocated
if (nuc->grid_.size() > 0) {
double max_E = nuc->grid_[0].energy.back();
int neutron = static_cast<int>(ParticleType::neutron);
int neutron = static_cast<int>(Particle::Type::neutron);
if (max_E == data::energy_max[neutron]) {
write_message("Maximum neutron transport energy: " +
std::to_string(data::energy_max[neutron]) + " eV for " +

9
src/dagmc.cpp
View file

@ -143,14 +143,14 @@ void load_dagmc_geometry()
c->universe_ = dagmc_univ_id; // set to zero for now
c->fill_ = C_NONE; // no fill, single universe
model::cells.push_back(c);
model::cells.emplace_back(c);
model::cell_map[c->id_] = i;
// Populate the Universe vector and dict
auto it = model::universe_map.find(dagmc_univ_id);
if (it == model::universe_map.end()) {
model::universes.push_back(new Universe());
model::universes.back()-> id_ = dagmc_univ_id;
model::universes.push_back(std::make_unique<Universe>());
model::universes.back()->id_ = dagmc_univ_id;
model::universes.back()->cells_.push_back(i);
model::universe_map[dagmc_univ_id] = model::universes.size() - 1;
} else {
@ -255,7 +255,6 @@ void load_dagmc_geometry()
// initialize surface objects
int n_surfaces = model::DAG->num_entities(2);
model::surfaces.resize(n_surfaces);
for (int i = 0; i < n_surfaces; i++) {
moab::EntityHandle surf_handle = model::DAG->entity_by_index(2, i+1);
@ -303,7 +302,7 @@ void load_dagmc_geometry()
}
// add to global array and map
model::surfaces[i] = s;
model::surfaces.emplace_back(s);
model::surface_map[s->id_] = s->id_;
}

4
src/eigenvalue.cpp
View file

@ -132,7 +132,7 @@ void synchronize_bank()
// Allocate temporary source bank
int64_t index_temp = 0;
std::vector<Bank> temp_sites(3*simulation::work);
std::vector<Particle::Bank> 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!");

233
src/geometry.cpp
View file

@ -34,21 +34,21 @@ std::vector<int64_t> overlap_check_count;
extern "C" bool
check_cell_overlap(Particle* p)
{
int n_coord = p->n_coord;
int n_coord = p->n_coord_;
// Loop through each coordinate level
for (int j = 0; j < n_coord; j++) {
Universe& univ = *model::universes[p->coord[j].universe];
Universe& univ = *model::universes[p->coord_[j].universe];
int n = univ.cells_.size();
// Loop through each cell on this level
for (auto index_cell : univ.cells_) {
Cell& c = *model::cells[index_cell];
if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) {
if (index_cell != p->coord[j].cell) {
if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) {
if (index_cell != p->coord_[j].cell) {
std::stringstream err_msg;
err_msg << "Overlapping cells detected: " << c.id_ << ", "
<< model::cells[p->coord[j].cell]->id_ << " on universe "
<< model::cells[p->coord_[j].cell]->id_ << " on universe "
<< univ.id_;
fatal_error(err_msg);
}
@ -74,36 +74,36 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
i_cell = *it;
// Make sure the search cell is in the same universe.
int i_universe = p->coord[p->n_coord-1].universe;
int i_universe = p->coord_[p->n_coord_-1].universe;
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->coord[p->n_coord-1].xyz};
Direction u {p->coord[p->n_coord-1].uvw};
auto surf = p->surface;
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord[p->n_coord-1].cell = i_cell;
p->coord_[p->n_coord_-1].cell = i_cell;
found = true;
break;
}
}
} else {
int i_universe = p->coord[p->n_coord-1].universe;
int i_universe = p->coord_[p->n_coord_-1].universe;
const auto& cells {model::universes[i_universe]->cells_};
for (auto it = cells.cbegin(); it != cells.cend(); it++) {
i_cell = *it;
// Make sure the search cell is in the same universe.
int i_universe = p->coord[p->n_coord-1].universe;
int i_universe = p->coord_[p->n_coord_-1].universe;
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->coord[p->n_coord-1].xyz};
Direction u {p->coord[p->n_coord-1].uvw};
auto surf = p->surface;
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord[p->n_coord-1].cell = i_cell;
p->coord_[p->n_coord_-1].cell = i_cell;
found = true;
break;
}
@ -126,37 +126,37 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
// Find the distribcell instance number.
if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) {
int offset = 0;
for (int i = 0; i < p->n_coord; i++) {
Cell& c_i {*model::cells[p->coord[i].cell]};
for (int i = 0; i < p->n_coord_; i++) {
const auto& c_i {*model::cells[p->coord_[i].cell]};
if (c_i.type_ == FILL_UNIVERSE) {
offset += c_i.offset_[c.distribcell_index_];
} else if (c_i.type_ == FILL_LATTICE) {
Lattice& lat {*model::lattices[p->coord[i+1].lattice-1]};
int i_xyz[3] {p->coord[i+1].lattice_x,
p->coord[i+1].lattice_y,
p->coord[i+1].lattice_z};
auto& lat {*model::lattices[p->coord_[i+1].lattice]};
int i_xyz[3] {p->coord_[i+1].lattice_x,
p->coord_[i+1].lattice_y,
p->coord_[i+1].lattice_z};
if (lat.are_valid_indices(i_xyz)) {
offset += lat.offset(c.distribcell_index_, i_xyz);
}
}
}
p->cell_instance = offset;
p->cell_instance_ = offset;
} else {
p->cell_instance = 0;
p->cell_instance_ = 0;
}
// Set the material and temperature.
p->last_material = p->material;
p->material_last_ = p->material_;
if (c.material_.size() > 1) {
p->material = c.material_[p->cell_instance];
p->material_ = c.material_[p->cell_instance_];
} else {
p->material = c.material_[0];
p->material_ = c.material_[0];
}
p->last_sqrtkT = p->sqrtkT;
p->sqrtkT_last_ = p->sqrtkT_;
if (c.sqrtkT_.size() > 1) {
p->sqrtkT = c.sqrtkT_[p->cell_instance];
p->sqrtkT_ = c.sqrtkT_[p->cell_instance_];
} else {
p->sqrtkT = c.sqrtkT_[0];
p->sqrtkT_ = c.sqrtkT_[0];
}
return true;
@ -166,44 +166,36 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
//! Found a lower universe, update this coord level then search the next.
// Set the lower coordinate level universe.
p->coord[p->n_coord].universe = c.fill_;
p->coord_[p->n_coord_].universe = c.fill_;
// Set the position and direction.
for (int i = 0; i < 3; i++) {
p->coord[p->n_coord].xyz[i] = p->coord[p->n_coord-1].xyz[i];
p->coord[p->n_coord].uvw[i] = p->coord[p->n_coord-1].uvw[i];
}
p->coord_[p->n_coord_].r = p->coord_[p->n_coord_-1].r;
p->coord_[p->n_coord_].u = p->coord_[p->n_coord_-1].u;
// Apply translation.
p->coord[p->n_coord].xyz[0] -= c.translation_.x;
p->coord[p->n_coord].xyz[1] -= c.translation_.y;
p->coord[p->n_coord].xyz[2] -= c.translation_.z;
p->coord_[p->n_coord_].r -= c.translation_;
// Apply rotation.
if (!c.rotation_.empty()) {
auto x = p->coord[p->n_coord].xyz[0];
auto y = p->coord[p->n_coord].xyz[1];
auto z = p->coord[p->n_coord].xyz[2];
p->coord[p->n_coord].xyz[0] = x*c.rotation_[3] + y*c.rotation_[4]
+ z*c.rotation_[5];
p->coord[p->n_coord].xyz[1] = x*c.rotation_[6] + y*c.rotation_[7]
+ z*c.rotation_[8];
p->coord[p->n_coord].xyz[2] = x*c.rotation_[9] + y*c.rotation_[10]
+ z*c.rotation_[11];
auto u = p->coord[p->n_coord].uvw[0];
auto v = p->coord[p->n_coord].uvw[1];
auto w = p->coord[p->n_coord].uvw[2];
p->coord[p->n_coord].uvw[0] = u*c.rotation_[3] + v*c.rotation_[4]
+ w*c.rotation_[5];
p->coord[p->n_coord].uvw[1] = u*c.rotation_[6] + v*c.rotation_[7]
+ w*c.rotation_[8];
p->coord[p->n_coord].uvw[2] = u*c.rotation_[9] + v*c.rotation_[10]
+ w*c.rotation_[11];
p->coord[p->n_coord].rotated = true;
Position r = p->coord_[p->n_coord_].r;
p->coord_[p->n_coord_].r.x = r.x*c.rotation_[3] + r.y*c.rotation_[4]
+ r.z*c.rotation_[5];
p->coord_[p->n_coord_].r.y = r.x*c.rotation_[6] + r.y*c.rotation_[7]
+ r.z*c.rotation_[8];
p->coord_[p->n_coord_].r.z = r.x*c.rotation_[9] + r.y*c.rotation_[10]
+ r.z*c.rotation_[11];
Direction u = p->coord_[p->n_coord_].u;
p->coord_[p->n_coord_].u.x = u.x*c.rotation_[3] + u.y*c.rotation_[4]
+ u.z*c.rotation_[5];
p->coord_[p->n_coord_].u.y = u.x*c.rotation_[6] + u.y*c.rotation_[7]
+ u.z*c.rotation_[8];
p->coord_[p->n_coord_].u.z = u.x*c.rotation_[9] + u.y*c.rotation_[10]
+ u.z*c.rotation_[11];
p->coord_[p->n_coord_].rotated = true;
}
// Update the coordinate level and recurse.
++p->n_coord;
++p->n_coord_;
return find_cell_inner(p, nullptr);
} else if (c.type_ == FILL_LATTICE) {
@ -213,35 +205,28 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
Lattice& lat {*model::lattices[c.fill_]};
// Determine lattice indices.
Position r {p->coord[p->n_coord-1].xyz};
Direction u {p->coord[p->n_coord-1].uvw};
r += TINY_BIT * u;
auto i_xyz = lat.get_indices(r);
auto i_xyz = lat.get_indices(p->r_local(), p->u_local());
// Store lower level coordinates.
r = lat.get_local_position(p->coord[p->n_coord-1].xyz, i_xyz);
p->coord[p->n_coord].xyz[0] = r.x;
p->coord[p->n_coord].xyz[1] = r.y;
p->coord[p->n_coord].xyz[2] = r.z;
p->coord[p->n_coord].uvw[0] = u.x;
p->coord[p->n_coord].uvw[1] = u.y;
p->coord[p->n_coord].uvw[2] = u.z;
Position r = lat.get_local_position(p->r_local(), i_xyz);
p->coord_[p->n_coord_].r = r;
p->coord_[p->n_coord_].u = p->u_local();
// Set lattice indices.
p->coord[p->n_coord].lattice = c.fill_ + 1;
p->coord[p->n_coord].lattice_x = i_xyz[0];
p->coord[p->n_coord].lattice_y = i_xyz[1];
p->coord[p->n_coord].lattice_z = i_xyz[2];
p->coord_[p->n_coord_].lattice = c.fill_;
p->coord_[p->n_coord_].lattice_x = i_xyz[0];
p->coord_[p->n_coord_].lattice_y = i_xyz[1];
p->coord_[p->n_coord_].lattice_z = i_xyz[2];
// Set the lower coordinate level universe.
if (lat.are_valid_indices(i_xyz)) {
p->coord[p->n_coord].universe = lat[i_xyz];
p->coord_[p->n_coord_].universe = lat[i_xyz];
} else {
if (lat.outer_ != NO_OUTER_UNIVERSE) {
p->coord[p->n_coord].universe = lat.outer_;
p->coord_[p->n_coord_].universe = lat.outer_;
} else {
std::stringstream err_msg;
err_msg << "Particle " << p->id << " is outside lattice "
err_msg << "Particle " << p->id_ << " is outside lattice "
<< lat.id_ << " but the lattice has no defined outer "
"universe.";
warning(err_msg);
@ -250,7 +235,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
}
// Update the coordinate level and recurse.
++p->n_coord;
++p->n_coord_;
return find_cell_inner(p, nullptr);
}
}
@ -264,22 +249,22 @@ extern "C" bool
find_cell(Particle* p, bool use_neighbor_lists)
{
// Determine universe (if not yet set, use root universe).
int i_universe = p->coord[p->n_coord-1].universe;
int i_universe = p->coord_[p->n_coord_-1].universe;
if (i_universe == C_NONE) {
p->coord[0].universe = model::root_universe;
p->n_coord = 1;
p->coord_[0].universe = model::root_universe;
p->n_coord_ = 1;
i_universe = model::root_universe;
}
// Reset all the deeper coordinate levels.
for (int i = p->n_coord; i < MAX_COORD; i++) {
p->coord[i].reset();
for (int i = p->n_coord_; i < MAX_COORD; i++) {
p->coord_[i].reset();
}
if (use_neighbor_lists) {
// Get the cell this particle was in previously.
auto coord_lvl = p->n_coord - 1;
auto i_cell = p->coord[coord_lvl].cell;
auto coord_lvl = p->n_coord_ - 1;
auto i_cell = p->coord_[coord_lvl].cell;
Cell& c {*model::cells[i_cell]};
// Search for the particle in that cell's neighbor list. Return if we
@ -291,7 +276,7 @@ find_cell(Particle* p, bool use_neighbor_lists)
// cells in this universe, and update the neighbor list if we find a new
// neighboring cell.
found = find_cell_inner(p, nullptr);
if (found) c.neighbors_.push_back(p->coord[coord_lvl].cell);
if (found) c.neighbors_.push_back(p->coord_[coord_lvl].cell);
return found;
} else {
@ -305,55 +290,52 @@ find_cell(Particle* p, bool use_neighbor_lists)
extern "C" void
cross_lattice(Particle* p, int lattice_translation[3])
{
Lattice& lat {*model::lattices[p->coord[p->n_coord-1].lattice-1]};
auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]};
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
msg << " Crossing lattice " << lat.id_ << ". Current position ("
<< p->coord[p->n_coord-1].lattice_x << ","
<< p->coord[p->n_coord-1].lattice_y << ","
<< p->coord[p->n_coord-1].lattice_z << ")";
<< p->coord_[p->n_coord_-1].lattice_x << ","
<< p->coord_[p->n_coord_-1].lattice_y << ","
<< p->coord_[p->n_coord_-1].lattice_z << ")";
write_message(msg, 1);
}
// Set the lattice indices.
p->coord[p->n_coord-1].lattice_x += lattice_translation[0];
p->coord[p->n_coord-1].lattice_y += lattice_translation[1];
p->coord[p->n_coord-1].lattice_z += lattice_translation[2];
std::array<int, 3> i_xyz {p->coord[p->n_coord-1].lattice_x,
p->coord[p->n_coord-1].lattice_y,
p->coord[p->n_coord-1].lattice_z};
p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0];
p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1];
p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2];
std::array<int, 3> i_xyz {p->coord_[p->n_coord_-1].lattice_x,
p->coord_[p->n_coord_-1].lattice_y,
p->coord_[p->n_coord_-1].lattice_z};
// Set the new coordinate position.
auto r = lat.get_local_position(p->coord[p->n_coord-2].xyz, i_xyz);
p->coord[p->n_coord-1].xyz[0] = r.x;
p->coord[p->n_coord-1].xyz[1] = r.y;
p->coord[p->n_coord-1].xyz[2] = r.z;
p->r_local() = lat.get_local_position(p->coord_[p->n_coord_-2].r, i_xyz);
if (!lat.are_valid_indices(i_xyz)) {
// The particle is outside the lattice. Search for it from the base coords.
p->n_coord = 1;
p->n_coord_ = 1;
bool found = find_cell(p, 0);
if (!found && p->alive) {
if (!found && p->alive_) {
std::stringstream err_msg;
err_msg << "Could not locate particle " << p->id
err_msg << "Could not locate particle " << p->id_
<< " after crossing a lattice boundary";
p->mark_as_lost(err_msg);
}
} else {
// Find cell in next lattice element.
p->coord[p->n_coord-1].universe = lat[i_xyz];
p->coord_[p->n_coord_-1].universe = lat[i_xyz];
bool found = find_cell(p, 0);
if (!found) {
// A particle crossing the corner of a lattice tile may not be found. In
// this case, search for it from the base coords.
p->n_coord = 1;
p->n_coord_ = 1;
bool found = find_cell(p, 0);
if (!found && p->alive) {
if (!found && p->alive_) {
std::stringstream err_msg;
err_msg << "Could not locate particle " << p->id
err_msg << "Could not locate particle " << p->id_
<< " after crossing a lattice boundary";
p->mark_as_lost(err_msg);
}
@ -377,21 +359,21 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
std::array<int, 3> level_lat_trans;
// Loop over each coordinate level.
for (int i = 0; i < p->n_coord; i++) {
Position r {p->coord[i].xyz};
Direction u {p->coord[i].uvw};
Cell& c {*model::cells[p->coord[i].cell]};
for (int i = 0; i < p->n_coord_; i++) {
Position r {p->coord_[i].r};
Direction u {p->coord_[i].u};
Cell& c {*model::cells[p->coord_[i].cell]};
// Find the oncoming surface in this cell and the distance to it.
auto surface_distance = c.distance(r, u, p->surface);
auto surface_distance = c.distance(r, u, p->surface_);
d_surf = surface_distance.first;
level_surf_cross = surface_distance.second;
// Find the distance to the next lattice tile crossing.
if (p->coord[i].lattice != F90_NONE) {
Lattice& lat {*model::lattices[p->coord[i].lattice-1]};
std::array<int, 3> i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y,
p->coord[i].lattice_z};
if (p->coord_[i].lattice != C_NONE) {
auto& lat {*model::lattices[p->coord_[i].lattice]};
std::array<int, 3> i_xyz {p->coord_[i].lattice_x, p->coord_[i].lattice_y,
p->coord_[i].lattice_z};
//TODO: refactor so both lattice use the same position argument (which
//also means the lat.type attribute can be removed)
std::pair<double, std::array<int, 3>> lattice_distance;
@ -400,8 +382,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
lattice_distance = lat.distance(r, u, i_xyz);
break;
case LatticeType::hex:
Position r_hex {p->coord[i-1].xyz[0], p->coord[i-1].xyz[1],
p->coord[i].xyz[2]};
Position r_hex {p->coord_[i-1].r.x, p->coord_[i-1].r.y,
p->coord_[i].r.z};
lattice_distance = lat.distance(r_hex, u, i_xyz);
break;
}
@ -410,7 +392,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
if (d_lat < 0) {
std::stringstream err_msg;
err_msg << "Particle " << p->id
err_msg << "Particle " << p->id_
<< " had a negative distance to a lattice boundary";
p->mark_as_lost(err_msg);
}
@ -466,23 +448,20 @@ extern "C" int
openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
{
Particle p;
p.initialize();
std::copy(xyz, xyz + 3, p.coord[0].xyz);
p.coord[0].uvw[0] = 0.0;
p.coord[0].uvw[1] = 0.0;
p.coord[0].uvw[2] = 1.0;
p.r() = Position{xyz};
p.u() = {0.0, 0.0, 1.0};
if (!find_cell(&p, false)) {
std::stringstream msg;
msg << "Could not find cell at position (" << xyz[0] << ", " << xyz[1]
<< ", " << xyz[2] << ").";
msg << "Could not find cell at position (" << p.r().x << ", " << p.r().y
<< ", " << p.r().z << ").";
set_errmsg(msg);
return OPENMC_E_GEOMETRY;
}
*index = p.coord[p.n_coord-1].cell;
*instance = p.cell_instance;
*index = p.coord_[p.n_coord_-1].cell;
*instance = p.cell_instance_;
return 0;
}

23
src/geometry_aux.cpp
View file

@ -66,7 +66,7 @@ void
adjust_indices()
{
// Adjust material/fill idices.
for (Cell* c : model::cells) {
for (auto& c : model::cells) {
if (c->fill_ != C_NONE) {
int32_t id = c->fill_;
auto search_univ = model::universe_map.find(id);
@ -102,7 +102,7 @@ adjust_indices()
}
// Change cell.universe values from IDs to indices.
for (Cell* c : model::cells) {
for (auto& c : model::cells) {
auto search = model::universe_map.find(c->universe_);
if (search != model::universe_map.end()) {
c->universe_ = search->second;
@ -115,7 +115,7 @@ adjust_indices()
}
// Change all lattice universe values from IDs to indices.
for (Lattice* l : model::lattices) {
for (auto& l : model::lattices) {
l->adjust_indices();
}
}
@ -125,7 +125,7 @@ adjust_indices()
void
assign_temperatures()
{
for (Cell* c : model::cells) {
for (auto& c : model::cells) {
// Ignore non-material cells and cells with defined temperature.
if (c->material_.size() == 0) continue;
if (c->sqrtkT_.size() > 0) continue;
@ -234,12 +234,12 @@ find_root_universe()
{
// Find all the universes listed as a cell fill.
std::unordered_set<int32_t> fill_univ_ids;
for (Cell* c : model::cells) {
for (const auto& c : model::cells) {
fill_univ_ids.insert(c->fill_);
}
// Find all the universes contained in a lattice.
for (Lattice* lat : model::lattices) {
for (const auto& lat : model::lattices) {
for (auto it = lat->begin(); it != lat->end(); ++it) {
fill_univ_ids.insert(*it);
}
@ -318,7 +318,7 @@ prepare_distribcell()
// unique distribcell array index.
int distribcell_index = 0;
std::vector<int32_t> target_univ_ids;
for (Universe* u : model::universes) {
for (const auto& u : model::universes) {
for (auto cell_indx : u->cells_) {
if (distribcells.find(cell_indx) != distribcells.end()) {
model::cells[cell_indx]->distribcell_index_ = distribcell_index;
@ -330,19 +330,19 @@ prepare_distribcell()
// Allocate the cell and lattice offset tables.
int n_maps = target_univ_ids.size();
for (Cell* c : model::cells) {
for (auto& c : model::cells) {
if (c->type_ != FILL_MATERIAL) {
c->offset_.resize(n_maps, C_NONE);
}
}
for (Lattice* lat : model::lattices) {
for (auto& lat : model::lattices) {
lat->allocate_offset_table(n_maps);
}
// Fill the cell and lattice offset tables.
for (int map = 0; map < target_univ_ids.size(); map++) {
auto target_univ_id = target_univ_ids[map];
for (Universe* univ : model::universes) {
for (const auto& univ : model::universes) {
int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation.
for (int32_t cell_indx : univ->cells_) {
Cell& c = *model::cells[cell_indx];
@ -525,15 +525,12 @@ maximum_levels(int32_t univ)
void
free_memory_geometry()
{
for (Cell* c : model::cells) {delete c;}
model::cells.clear();
model::cell_map.clear();
for (Universe* u : model::universes) {delete u;}
model::universes.clear();
model::universe_map.clear();
for (Lattice* lat : model::lattices) {delete lat;}
model::lattices.clear();
model::lattice_map.clear();

12
src/initialize.cpp
View file

@ -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);

69
src/lattice.cpp
View file

@ -19,10 +19,8 @@ namespace openmc {
//==============================================================================
namespace model {
std::vector<Lattice*> lattices;
std::unordered_map<int32_t, int32_t> lattice_map;
std::vector<std::unique_ptr<Lattice>> lattices;
std::unordered_map<int32_t, int32_t> lattice_map;
}
//==============================================================================
@ -287,15 +285,38 @@ const
//==============================================================================
std::array<int, 3>
RectLattice::get_indices(Position r) const
RectLattice::get_indices(Position r, Direction u) const
{
int ix {static_cast<int>(std::ceil((r.x - lower_left_.x) / pitch_.x))-1};
int iy {static_cast<int>(std::ceil((r.y - lower_left_.y) / pitch_.y))-1};
int iz;
if (is_3d_) {
iz = static_cast<int>(std::ceil((r.z - lower_left_.z) / pitch_.z))-1;
// Determine x index, accounting for coincidence
double ix_ {(r.x - lower_left_.x) / pitch_.x};
long ix_close {std::lround(ix_)};
int ix;
if (std::abs(ix_ - ix_close) < FP_COINCIDENT) {
ix = (u.x > 0) ? ix_close : ix_close - 1;
} else {
iz = 0;
ix = std::floor(ix_);
}
// Determine y index, accounting for coincidence
double iy_ {(r.y - lower_left_.y) / pitch_.y};
long iy_close {std::lround(iy_)};
int iy;
if (std::abs(iy_ - iy_close) < FP_COINCIDENT) {
iy = (u.y > 0) ? iy_close : iy_close - 1;
} else {
iy = std::floor(iy_);
}
// Determine z index, accounting for coincidence
int iz = 0;
if (is_3d_) {
double iz_ {(r.z - lower_left_.z) / pitch_.z};
long iz_close {std::lround(iz_)};
if (std::abs(iz_ - iz_close) < FP_COINCIDENT) {
iz = (u.z > 0) ? iz_close : iz_close - 1;
} else {
iz = std::floor(iz_);
}
}
return {ix, iy, iz};
}
@ -687,8 +708,13 @@ const
//==============================================================================
std::array<int, 3>
HexLattice::get_indices(Position r) const
HexLattice::get_indices(Position r, Direction u) const
{
// The implementation for HexLattice currently doesn't use direction
// information. As a result, we move the position slightly forward to
// determine what lattice index the particle is most likely to be in.
r += TINY_BIT * u;
// Offset the xyz by the lattice center.
Position r_o {r.x - center_.x, r.y - center_.y, r.z};
if (is_3d_) {r_o.z -= center_.z;}
@ -696,7 +722,7 @@ HexLattice::get_indices(Position r) const
// Index the z direction.
std::array<int, 3> out;
if (is_3d_) {
out[2] = static_cast<int>(std::ceil(r_o.z / pitch_[1] + 0.5 * n_axial_))-1;
out[2] = std::floor(r_o.z / pitch_[1] + 0.5 * n_axial_);
} else {
out[2] = 0;
}
@ -704,9 +730,8 @@ HexLattice::get_indices(Position r) const
// Convert coordinates into skewed bases. The (x, alpha) basis is used to
// find the index of the global coordinates to within 4 cells.
double alpha = r_o.y - r_o.x / std::sqrt(3.0);
out[0] = static_cast<int>(std::floor(r_o.x
/ (0.5*std::sqrt(3.0) * pitch_[0])));
out[1] = static_cast<int>(std::floor(alpha / pitch_[0]));
out[0] = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
out[1] = std::floor(alpha / pitch_[0]);
// Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
// the array is offset so that the indices never go below 0).
@ -739,12 +764,12 @@ HexLattice::get_indices(Position r) const
// Select the minimum squared distance which corresponds to the cell the
// coordinates are in.
if (k_min == 2) {
out[0] += 1;
++out[0];
} else if (k_min == 3) {
out[1] += 1;
++out[1];
} else if (k_min == 4) {
out[0] += 1;
out[1] += 1;
++out[0];
++out[1];
}
return out;
@ -867,10 +892,10 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
void read_lattices(pugi::xml_node node)
{
for (pugi::xml_node lat_node : node.children("lattice")) {
model::lattices.push_back(new RectLattice(lat_node));
model::lattices.push_back(std::make_unique<RectLattice>(lat_node));
}
for (pugi::xml_node lat_node : node.children("hex_lattice")) {
model::lattices.push_back(new HexLattice(lat_node));
model::lattices.push_back(std::make_unique<HexLattice>(lat_node));
}
// Fill the lattice map.

36
src/material.cpp
View file

@ -37,7 +37,7 @@ namespace openmc {
namespace model {
std::vector<Material*> materials;
std::vector<std::unique_ptr<Material>> materials;
std::unordered_map<int32_t, int32_t> material_map;
} // namespace model
@ -652,19 +652,18 @@ void Material::calculate_xs(const Particle& p) const
simulation::material_xs.fission = 0.0;
simulation::material_xs.nu_fission = 0.0;
if (p.type == static_cast<int>(ParticleType::neutron)) {
if (p.type_ == Particle::Type::neutron) {
this->calculate_neutron_xs(p);
} else if (p.type == static_cast<int>(ParticleType::photon)) {
} else if (p.type_ == Particle::Type::photon) {
this->calculate_photon_xs(p);
}
}
void Material::calculate_neutron_xs(const Particle& p) const
{
int neutron = static_cast<int>(ParticleType::neutron);
// Find energy index on energy grid
int i_grid = std::log(p.E/data::energy_min[neutron])/simulation::log_spacing;
int neutron = static_cast<int>(Particle::Type::neutron);
int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing;
// Determine if this material has S(a,b) tables
bool check_sab = (thermal_tables_.size() > 0);
@ -691,7 +690,7 @@ void Material::calculate_neutron_xs(const Particle& p) const
// If particle energy is greater than the highest energy for the
// S(a,b) table, then don't use the S(a,b) table
if (p.E > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
if (p.E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
// Increment position in thermal_tables_
++j;
@ -709,12 +708,12 @@ void Material::calculate_neutron_xs(const Particle& p) const
// Calculate microscopic cross section for this nuclide
const auto& micro {simulation::micro_xs[i_nuclide]};
if (p.E != micro.last_E
|| p.sqrtkT != micro.last_sqrtkT
if (p.E_ != micro.last_E
|| p.sqrtkT_ != micro.last_sqrtkT
|| i_sab != micro.index_sab
|| sab_frac != micro.sab_frac) {
data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E, i_grid,
p.sqrtkT, sab_frac);
data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E_, i_grid,
p.sqrtkT_, sab_frac);
}
// ======================================================================
@ -748,8 +747,8 @@ void Material::calculate_photon_xs(const Particle& p) const
// Calculate microscopic cross section for this nuclide
const auto& micro {simulation::micro_photon_xs[i_element]};
if (p.E != micro.last_E) {
data::elements[i_element].calculate_xs(p.E);
if (p.E_ != micro.last_E) {
data::elements[i_element].calculate_xs(p.E_);
}
// ========================================================================
@ -895,7 +894,7 @@ void read_materials_xml()
// Loop over XML material elements and populate the array.
pugi::xml_node root = doc.document_element();
for (pugi::xml_node material_node : root.children("material")) {
model::materials.push_back(new Material(material_node));
model::materials.push_back(std::make_unique<Material>(material_node));
}
model::materials.shrink_to_fit();
@ -915,7 +914,6 @@ void read_materials_xml()
void free_memory_material()
{
for (Material *mat : model::materials) {delete mat;}
model::materials.clear();
model::material_map.clear();
}
@ -942,7 +940,7 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density)
{
int err = 0;
if (index >= 0 && index < model::materials.size()) {
Material* m = model::materials[index];
auto& m = model::materials[index];
// Check if nuclide is already in material
for (int i = 0; i < m->nuclide_.size(); ++i) {
@ -1032,7 +1030,7 @@ extern "C" int
openmc_material_get_volume(int32_t index, double* volume)
{
if (index >= 0 && index < model::materials.size()) {
Material* m = model::materials[index];
auto& m = model::materials[index];
if (m->volume_ >= 0.0) {
*volume = m->volume_;
return 0;
@ -1112,7 +1110,7 @@ extern "C" int
openmc_material_set_volume(int32_t index, double volume)
{
if (index >= 0 && index < model::materials.size()) {
Material* m = model::materials[index];
auto& m {model::materials[index]};
if (volume >= 0.0) {
m->volume_ = volume;
return 0;
@ -1132,7 +1130,7 @@ openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end)
if (index_start) *index_start = model::materials.size();
if (index_end) *index_end = model::materials.size() + n - 1;
for (int32_t i = 0; i < n; i++) {
model::materials.push_back(new Material());
model::materials.push_back(std::make_unique<Material>());
}
return 0;
}

68
src/math_functions.cpp
View file

@ -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();

20
src/mesh.cpp
View file

@ -385,9 +385,9 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
// just a bit for the purposes of determining if there was an intersection
// in case the mesh surfaces coincide with lattice/geometric surfaces which
// might produce finite-precision errors.
Position last_r {p->last_xyz};
Position r {p->coord[0].xyz};
Direction u {p->coord[0].uvw};
Position last_r {p->r_last_};
Position r {p->r()};
Direction u {p->u()};
Position r0 = last_r + TINY_BIT*u;
Position r1 = r - TINY_BIT*u;
@ -542,9 +542,9 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins
// Determine if the track intersects the tally mesh.
// Copy the starting and ending coordinates of the particle.
Position r0 {p->last_xyz_current};
Position r1 {p->coord[0].xyz};
Direction u {p->coord[0].uvw};
Position r0 {p->r_last_current_};
Position r1 {p->r()};
Direction u {p->u()};
// Determine indices for starting and ending location.
int n = n_dimension_;
@ -695,7 +695,7 @@ void RegularMesh::to_hdf5(hid_t group) const
close_group(mesh_group);
}
xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* bank,
int n_energy, const double* energies, bool* outside) const
{
// Determine shape of array for counts
@ -713,7 +713,7 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
for (int64_t i = 0; i < n; ++i) {
// determine scoring bin for entropy mesh
int mesh_bin = get_bin({bank[i].xyz});
int mesh_bin = get_bin(bank[i].r);
// if outside mesh, skip particle
if (mesh_bin < 0) {
@ -771,7 +771,7 @@ openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
{
if (index_start) *index_start = model::meshes.size();
for (int i = 0; i < n; ++i) {
model::meshes.emplace_back(new RegularMesh{});
model::meshes.push_back(std::make_unique<RegularMesh>());
}
if (index_end) *index_end = model::meshes.size() - 1;
@ -909,7 +909,7 @@ void read_meshes(pugi::xml_node root)
{
for (auto node : root.children("mesh")) {
// Read mesh and add to vector
model::meshes.emplace_back(new RegularMesh{node});
model::meshes.push_back(std::make_unique<RegularMesh>(node));
// Map ID to position in vector
model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;

22
src/mgxs.cpp
View file

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

8
src/mgxs_interface.cpp
View file

@ -136,7 +136,7 @@ void create_macro_xs()
for (int i = 0; i < model::materials.size(); ++i) {
if (kTs[i].size() > 0) {
// Convert atom_densities to a vector
Material* mat = model::materials[i];
auto& mat {model::materials[i]};
std::vector<double> atom_densities(mat->atom_density_.begin(),
mat->atom_density_.end());
@ -235,7 +235,7 @@ void read_mg_cross_sections_header()
}
// Get the minimum and maximum energies
int neutron = static_cast<int>(ParticleType::neutron);
int neutron = static_cast<int>(Particle::Type::neutron);
data::energy_min[neutron] = data::energy_bins.back();
data::energy_max[neutron] = data::energy_bins.front();
@ -248,10 +248,10 @@ void read_mg_cross_sections_header()
//==============================================================================
void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
calculate_xs_c(int i_mat, int gin, double sqrtkT, Direction u,
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, u, total_xs, abs_xs,
nu_fiss_xs);
}

6
src/nuclide.cpp
View file

@ -287,7 +287,7 @@ void Nuclide::create_derived()
auto xs = xt::adapt(rx->xs_[t].value);
for (const auto& p : rx->products_) {
if (p.particle_ == ParticleType::photon) {
if (p.particle_ == Particle::Type::photon) {
auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD);
for (int k = 0; k < n; ++k) {
double E = grid_[t].energy[k+j];
@ -391,7 +391,7 @@ void Nuclide::create_derived()
void Nuclide::init_grid()
{
int neutron = static_cast<int>(ParticleType::neutron);
int neutron = static_cast<int>(Particle::Type::neutron);
double E_min = data::energy_min[neutron];
double E_max = data::energy_max[neutron];
int M = settings::n_log_bins;
@ -440,7 +440,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
for (int i = 1; i < rx->products_.size(); ++i) {
// Skip any non-neutron products
const auto& product = rx->products_[i];
if (product.particle_ != ParticleType::neutron) continue;
if (product.particle_ != Particle::Type::neutron) continue;
// Evaluate yield
if (product.emission_mode_ == EmissionMode::delayed) {

54
src/output.cpp
View file

@ -142,64 +142,64 @@ std::string time_stamp()
extern "C" void print_particle(Particle* p)
{
// Display particle type and ID.
switch (p->type) {
case static_cast<int>(ParticleType::neutron):
switch (p->type_) {
case Particle::Type::neutron:
std::cout << "Neutron ";
break;
case static_cast<int>(ParticleType::photon):
case Particle::Type::photon:
std::cout << "Photon ";
break;
case static_cast<int>(ParticleType::electron):
case Particle::Type::electron:
std::cout << "Electron ";
break;
case static_cast<int>(ParticleType::positron):
case Particle::Type::positron:
std::cout << "Positron ";
break;
default:
std::cout << "Unknown Particle ";
}
std::cout << p->id << "\n";
std::cout << p->id_ << "\n";
// Display particle geometry hierarchy.
for (auto i = 0; i < p->n_coord; i++) {
for (auto i = 0; i < p->n_coord_; i++) {
std::cout << " Level " << i << "\n";
if (p->coord[i].cell != C_NONE) {
const Cell& c {*model::cells[p->coord[i].cell]};
if (p->coord_[i].cell != C_NONE) {
const Cell& c {*model::cells[p->coord_[i].cell]};
std::cout << " Cell = " << c.id_ << "\n";
}
if (p->coord[i].universe != C_NONE) {
const Universe& u {*model::universes[p->coord[i].universe]};
if (p->coord_[i].universe != C_NONE) {
const Universe& u {*model::universes[p->coord_[i].universe]};
std::cout << " Universe = " << u.id_ << "\n";
}
if (p->coord[i].lattice != F90_NONE) {
const Lattice& lat {*model::lattices[p->coord[i].lattice]};
if (p->coord_[i].lattice != C_NONE) {
const Lattice& lat {*model::lattices[p->coord_[i].lattice]};
std::cout << " Lattice = " << lat.id_ << "\n";
std::cout << " Lattice position = (" << p->coord[i].lattice_x
<< "," << p->coord[i].lattice_y << ","
<< p->coord[i].lattice_z << ")\n";
std::cout << " Lattice position = (" << p->coord_[i].lattice_x
<< "," << p->coord_[i].lattice_y << ","
<< p->coord_[i].lattice_z << ")\n";
}
std::cout << " xyz = " << p->coord[i].xyz[0] << " "
<< p->coord[i].xyz[1] << " " << p->coord[i].xyz[2] << "\n";
std::cout << " uvw = " << p->coord[i].uvw[0] << " "
<< p->coord[i].uvw[1] << " " << p->coord[i].uvw[2] << "\n";
std::cout << " r = (" << p->coord_[i].r.x << ", "
<< p->coord_[i].r.y << ", " << p->coord_[i].r.z << ")\n";
std::cout << " u = (" << p->coord_[i].u.x << ", "
<< p->coord_[i].u.y << ", " << p->coord_[i].u.z << ")\n";
}
// Display miscellaneous info.
if (p->surface != ERROR_INT) {
const Surface& surf {*model::surfaces[std::abs(p->surface)-1]};
std::cout << " Surface = " << std::copysign(surf.id_, p->surface) << "\n";
if (p->surface_ != ERROR_INT) {
const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]};
std::cout << " Surface = " << std::copysign(surf.id_, p->surface_) << "\n";
}
std::cout << " Weight = " << p->wgt << "\n";
std::cout << " Weight = " << p->wgt_ << "\n";
if (settings::run_CE) {
std::cout << " Energy = " << p->E << "\n";
std::cout << " Energy = " << p->E_ << "\n";
} else {
std::cout << " Energy Group = " << p->g << "\n";
std::cout << " Energy Group = " << p->g_ << "\n";
}
std::cout << " Delayed Group = " << p->delayed_group << "\n";
std::cout << " Delayed Group = " << p->delayed_group_ << "\n";
std::cout << "\n";
}

350
src/particle.cpp
View file

@ -37,7 +37,7 @@ LocalCoord::reset()
{
cell = C_NONE;
universe = C_NONE;
lattice = 0;
lattice = C_NONE;
lattice_x = 0;
lattice_y = 0;
rotated = false;
@ -47,90 +47,70 @@ LocalCoord::reset()
// Particle implementation
//==============================================================================
void
Particle::clear()
{
// reset any coordinate levels
for (int i=0; i<MAX_COORD; ++i) coord[i].reset();
}
void
Particle::create_secondary(const double* uvw, double E, int type, bool run_CE)
{
if (n_secondary == MAX_SECONDARY) {
fatal_error("Too many secondary particles created.");
}
int64_t n = n_secondary;
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;
n_secondary += 1;
}
void
Particle::initialize()
Particle::Particle()
{
// Clear coordinate lists
clear();
// Set particle to neutron that's alive
type = static_cast<int>(ParticleType::neutron);
alive = true;
// clear attributes
surface = 0;
cell_born = C_NONE;
material = C_NONE;
last_material = C_NONE;
last_sqrtkT = 0;
wgt = 1.0;
last_wgt = 1.0;
absorb_wgt = 0.0;
n_bank = 0;
wgt_bank = 0.0;
sqrtkT = -1.0;
n_collision = 0;
fission = false;
delayed_group = 0;
for (int i=0; i<MAX_DELAYED_GROUPS; ++i) {
n_delayed_bank[i] = 0;
for (int& n : n_delayed_bank_) {
n = 0;
}
g = 0;
}
// Set up base level coordinates
coord[0].universe = C_NONE;
n_coord = 1;
last_n_coord = 1;
void
Particle::clear()
{
// reset any coordinate levels
for (auto& level : coord_) level.reset();
n_coord_ = 1;
}
void
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 = type;
secondary_bank_[n].wgt = wgt_;
secondary_bank_[n].r = this->r();
secondary_bank_[n].u = u;
secondary_bank_[n].E = settings::run_CE ? E : g_;
++n_secondary_;
}
void
Particle::from_source(const Bank* src)
{
// set defaults
initialize();
// reset some attributes
this->clear();
alive_ = true;
surface_ = 0;
cell_born_ = C_NONE;
material_ = C_NONE;
n_collision_ = 0;
fission_ = false;
// copy attributes from source bank site
type = src->particle;
wgt = src->wgt;
last_wgt = src->wgt;
std::copy(src->xyz, src->xyz + 3, coord[0].xyz);
std::copy(src->uvw, src->uvw + 3, coord[0].uvw);
std::copy(src->xyz, src->xyz + 3, last_xyz_current);
std::copy(src->xyz, src->xyz + 3, last_xyz);
std::copy(src->uvw, src->uvw + 3, last_uvw);
type_ = src->particle;
wgt_ = src->wgt;
wgt_last_ = src->wgt;
this->r() = src->r;
this->u() = src->u;
r_last_current_ = src->r;
r_last_ = src->r;
u_last_ = src->u;
if (settings::run_CE) {
E = src->E;
g = 0;
E_ = src->E;
g_ = 0;
} else {
g = static_cast<int>(src->E);
last_g = static_cast<int>(src->E);
E = data::energy_bin_avg[g - 1];
g_ = static_cast<int>(src->E);
g_last_ = static_cast<int>(src->E);
E_ = data::energy_bin_avg[g_ - 1];
}
last_E = E;
E_last_ = E_;
}
void
@ -138,7 +118,7 @@ Particle::transport()
{
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || simulation::trace) {
write_message("Simulating Particle " + std::to_string(id));
write_message("Simulating Particle " + std::to_string(id_));
}
// Initialize number of events to zero
@ -146,7 +126,7 @@ Particle::transport()
// Add paricle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += wgt;
simulation::total_weight += wgt_;
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
@ -156,62 +136,62 @@ Particle::transport()
}
// Prepare to write out particle track.
if (write_track) add_particle_track();
if (write_track_) add_particle_track();
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty()) zero_flux_derivs();
while (true) {
// Set the random number stream
if (type == static_cast<int>(ParticleType::neutron)) {
if (type_ == Particle::Type::neutron) {
prn_set_stream(STREAM_TRACKING);
} else {
prn_set_stream(STREAM_PHOTON);
}
// Store pre-collision particle properties
last_wgt = wgt;
last_E = E;
std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw);
std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz);
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
// If the cell hasn't been determined based on the particle's location,
// initiate a search for the current cell. This generally happens at the
// beginning of the history and again for any secondary particles
if (coord[n_coord - 1].cell == C_NONE) {
if (coord_[n_coord_ - 1].cell == C_NONE) {
if (!find_cell(this, false)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id));
+ std::to_string(id_));
return;
}
// set birth cell attribute
if (cell_born == C_NONE) cell_born = coord[n_coord - 1].cell;
if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
}
// Write particle track.
if (write_track) write_particle_track(*this);
if (write_track_) write_particle_track(*this);
if (settings::check_overlaps) check_cell_overlap(this);
// Calculate microscopic and macroscopic cross sections
if (material != MATERIAL_VOID) {
if (material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (material != last_material || sqrtkT != last_sqrtkT) {
if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
// If the material is the same as the last material and the
// temperature hasn't changed, we don't need to lookup cross
// sections again.
model::materials[material]->calculate_xs(*this);
model::materials[material_]->calculate_xs(*this);
}
} else {
// Get the MG data
calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw,
calculate_xs_c(material_, g_, sqrtkT_, this->u_local(),
simulation::material_xs.total, simulation::material_xs.absorption,
simulation::material_xs.nu_fission);
// Finally, update the particle group while we have already checked
// for if multi-group
last_g = g;
g_last_ = g_;
}
} else {
simulation::material_xs.total = 0.0;
@ -230,8 +210,8 @@ Particle::transport()
// Sample a distance to collision
double d_collision;
if (type == static_cast<int>(ParticleType::electron) ||
type == static_cast<int>(ParticleType::positron)) {
if (type_ == Particle::Type::electron ||
type_ == Particle::Type::positron) {
d_collision = 0.0;
} else if (simulation::material_xs.total == 0.0) {
d_collision = INFINITY;
@ -243,11 +223,8 @@ Particle::transport()
double distance = std::min(d_boundary, d_collision);
// Advance particle
for (int j = 0; j < n_coord; ++j) {
// TODO: use Position
coord[j].xyz[0] += distance * coord[j].uvw[0];
coord[j].xyz[1] += distance * coord[j].uvw[1];
coord[j].xyz[2] += distance * coord[j].uvw[2];
for (int j = 0; j < n_coord_; ++j) {
coord_[j].r += distance * coord_[j].u;
}
// Score track-length tallies
@ -257,8 +234,8 @@ Particle::transport()
// Score track-length estimate of k-eff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type == static_cast<int>(ParticleType::neutron)) {
global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission;
type_ == Particle::Type::neutron) {
global_tally_tracklength += wgt_ * distance * simulation::material_xs.nu_fission;
}
// Score flux derivative accumulators for differential tallies.
@ -270,25 +247,25 @@ Particle::transport()
// ====================================================================
// PARTICLE CROSSES SURFACE
if (next_level > 0) n_coord = next_level;
if (next_level > 0) n_coord_ = next_level;
// Saving previous cell data
for (int j = 0; j < n_coord; ++j) {
last_cell[j] = coord[j].cell;
for (int j = 0; j < n_coord_; ++j) {
cell_last_[j] = coord_[j].cell;
}
last_n_coord = n_coord;
n_coord_last_ = n_coord_;
if (lattice_translation[0] != 0 || lattice_translation[1] != 0 ||
lattice_translation[2] != 0) {
// Particle crosses lattice boundary
surface = ERROR_INT;
surface_ = ERROR_INT;
cross_lattice(this, lattice_translation);
event = EVENT_LATTICE;
event_ = EVENT_LATTICE;
} else {
// Particle crosses surface
surface = surface_crossed;
surface_ = surface_crossed;
this->cross_surface();
event = EVENT_SURFACE;
event_ = EVENT_SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
@ -300,8 +277,8 @@ Particle::transport()
// Score collision estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type == static_cast<int>(ParticleType::neutron)) {
global_tally_collision += wgt * simulation::material_xs.nu_fission
type_ == Particle::Type::neutron) {
global_tally_collision += wgt_ * simulation::material_xs.nu_fission
/ simulation::material_xs.total;
}
@ -313,7 +290,7 @@ Particle::transport()
score_surface_tally(this, model::active_meshsurf_tallies);
// Clear surface component
surface = ERROR_INT;
surface_ = ERROR_INT;
if (settings::run_CE) {
collision(this);
@ -334,33 +311,33 @@ Particle::transport()
}
// Reset banked weight during collision
n_bank = 0;
wgt_bank = 0.0;
for (int& v : n_delayed_bank) v = 0;
n_bank_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
// Reset fission logical
fission = false;
fission_ = false;
// Save coordinates for tallying purposes
std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current);
r_last_current_ = this->r();
// Set last material to none since cross sections will need to be
// re-evaluated
last_material = C_NONE;
material_last_ = C_NONE;
// Set all uvws to base level -- right now, after a collision, only the
// base level uvws are changed
for (int j = 0; j < n_coord - 1; ++j) {
if (coord[j + 1].rotated) {
// Set all directions to base level -- right now, after a collision, only
// the base level directions are changed
for (int j = 0; j < n_coord_ - 1; ++j) {
if (coord_[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord[j].cell]->rotation_};
Direction u {coord[j].uvw};
coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z;
const auto& m {model::cells[coord_[j].cell]->rotation_};
const auto& u {coord_[j].u};
coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
} else {
// Otherwise, copy this level's direction
std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw);
coord_[j+1].u = coord_[j].u;
}
}
@ -371,27 +348,27 @@ Particle::transport()
// If particle has too many events, display warning and kill it
++n_event;
if (n_event == MAX_EVENTS) {
warning("Particle " + std::to_string(id) +
warning("Particle " + std::to_string(id_) +
" underwent maximum number of events.");
alive = false;
alive_ = false;
}
// Check for secondary particles if this particle is dead
if (!alive) {
if (!alive_) {
// If no secondary particles, break out of event loop
if (n_secondary == 0) break;
if (n_secondary_ == 0) break;
this->from_source(&secondary_bank[n_secondary - 1]);
--n_secondary;
this->from_source(&secondary_bank_[n_secondary_ - 1]);
--n_secondary_;
n_event = 0;
// Enter new particle in particle track file
if (write_track) add_particle_track();
if (write_track_) add_particle_track();
}
}
// Finish particle track output.
if (write_track) {
if (write_track_) {
write_particle_track(*this);
finalize_particle_track(*this);
}
@ -400,9 +377,9 @@ Particle::transport()
void
Particle::cross_surface()
{
int i_surface = std::abs(surface);
int i_surface = std::abs(surface_);
// TODO: off-by-one
const auto& surf {model::surfaces[i_surface - 1]};
const auto& surf {model::surfaces[i_surface - 1].get()};
if (settings::verbosity >= 10 || simulation::trace) {
write_message(" Crossing surface " + std::to_string(surf->id_));
}
@ -412,7 +389,7 @@ Particle::cross_surface()
// PARTICLE LEAKS OUT OF PROBLEM
// Kill particle
alive = false;
alive_ = false;
// Score any surface current tallies -- note that the particle is moved
// forward slightly so that if the mesh boundary is on the surface, it is
@ -422,15 +399,12 @@ Particle::cross_surface()
// TODO: Find a better solution to score surface currents than
// physically moving the particle forward slightly
// TODO: Use Position
coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0];
coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1];
coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2];
this->r() += TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
}
// Score to global leakage tally
global_tally_leakage += wgt;
global_tally_leakage += wgt_;
// Display message
if (settings::verbosity >= 10 || simulation::trace) {
@ -443,8 +417,8 @@ Particle::cross_surface()
// PARTICLE REFLECTS FROM SURFACE
// Do not handle reflective boundary conditions on lower universes
if (n_coord != 1) {
this->mark_as_lost("Cannot reflect particle " + std::to_string(id) +
if (n_coord_ != 1) {
this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) +
" off surface in a lower universe.");
return;
}
@ -462,32 +436,27 @@ Particle::cross_surface()
if (!model::active_meshsurf_tallies.empty()) {
Position r {coord[0].xyz};
coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0];
coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1];
coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2];
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
std::copy(&r.x, &r.x + 3, coord[0].xyz);
this->r() = r;
}
// Reflect particle off surface
Direction u = surf->reflect(coord[0].xyz, coord[0].uvw);
Direction u = surf->reflect(this->r(), this->u());
// Make sure new particle direction is normalized
double norm = u.norm();
coord[0].uvw[0] = u.x/norm;
coord[0].uvw[1] = u.y/norm;
coord[0].uvw[2] = u.z/norm;
this->u() = u / u.norm();
// Reassign particle's cell and surface
coord[0].cell = last_cell[last_n_coord - 1];
surface = -surface;
coord_[0].cell = cell_last_[n_coord_last_ - 1];
surface_ = -surface_;
// If a reflective surface is coincident with a lattice or universe
// boundary, it is necessary to redetermine the particle's coordinates in
// the lower universes.
n_coord = 1;
n_coord_ = 1;
if (!find_cell(this, true)) {
this->mark_as_lost("Couldn't find particle after reflecting from surface "
+ std::to_string(surf->id_) + ".");
@ -495,9 +464,7 @@ Particle::cross_surface()
}
// Set previous coordinate going slightly past surface crossing
last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT*coord[0].uvw[0];
last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT*coord[0].uvw[1];
last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT*coord[0].uvw[2];
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
@ -510,8 +477,8 @@ Particle::cross_surface()
// PERIODIC BOUNDARY
// Do not handle periodic boundary conditions on lower universes
if (n_coord != 1) {
this->mark_as_lost("Cannot transfer particle " + std::to_string(id) +
if (n_coord_ != 1) {
this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) +
" across surface in a lower universe. Boundary conditions must be "
"applied to root universe.");
return;
@ -521,34 +488,28 @@ Particle::cross_surface()
// particle to change -- artificially move the particle slightly back in
// case the surface crossing is coincident with a mesh boundary
if (!model::active_meshsurf_tallies.empty()) {
Position r {coord[0].xyz};
coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0];
coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1];
coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2];
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
std::copy(&r.x, &r.x + 3, coord[0].xyz);
this->r() = r;
}
// Get a pointer to the partner periodic surface
auto surf_p = dynamic_cast<PeriodicSurface*>(surf);
auto other = dynamic_cast<PeriodicSurface*>(
model::surfaces[surf_p->i_periodic_]);
model::surfaces[surf_p->i_periodic_].get());
// Adjust the particle's location and direction.
Position r {coord[0].xyz};
Direction u {coord[0].uvw};
bool rotational = other->periodic_translate(surf_p, r, u);
std::copy(&r.x, &r.x + 3, coord[0].xyz);
std::copy(&u.x, &u.x + 3, coord[0].uvw);
bool rotational = other->periodic_translate(surf_p, this->r(), this->u());
// Reassign particle's surface
// TODO: off-by-one
surface = rotational ?
surface_ = rotational ?
surf_p->i_periodic_ + 1 :
std::copysign(surf_p->i_periodic_ + 1, surface);
std::copysign(surf_p->i_periodic_ + 1, surface_);
// Figure out what cell particle is in now
n_coord = 1;
n_coord_ = 1;
if (!find_cell(this, true)) {
this->mark_as_lost("Couldn't find particle after hitting periodic "
@ -557,9 +518,7 @@ Particle::cross_surface()
}
// Set previous coordinate going slightly past surface crossing
last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT * coord[0].uvw[0];
last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT * coord[0].uvw[1];
last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT * coord[0].uvw[2];
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
@ -574,18 +533,18 @@ Particle::cross_surface()
#ifdef DAGMC
if (settings::dagmc) {
auto cellp = dynamic_cast<DAGCell*>(model::cells[last_cell[0]]);
auto cellp = dynamic_cast<DAGCell*>(model::cells[cell_last_[0]].get());
// TODO: off-by-one
auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface) - 1]);
auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface_) - 1].get());
int32_t i_cell = next_cell(cellp, surfp) - 1;
// save material and temp
last_material = material;
last_sqrtkT = sqrtkT;
material_last_ = material_;
sqrtkT_last_ = sqrtkT_;
// set new cell value
coord[0].cell = i_cell;
cell_instance = 0;
material = model::cells[i_cell]->material_[0];
sqrtkT = model::cells[i_cell]->sqrtkT_[0];
coord_[0].cell = i_cell;
cell_instance_ = 0;
material_ = model::cells[i_cell]->material_[0];
sqrtkT_ = model::cells[i_cell]->sqrtkT_[0];
return;
}
#endif
@ -596,8 +555,8 @@ Particle::cross_surface()
// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
// Remove lower coordinate levels and assignment of surface
surface = ERROR_INT;
n_coord = 1;
surface_ = ERROR_INT;
n_coord_ = 1;
bool found = find_cell(this, false);
if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) {
@ -606,16 +565,14 @@ Particle::cross_surface()
// the particle is really traveling tangent to a surface, if we move it
// forward a tiny bit it should fix the problem.
n_coord = 1;
coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0];
coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1];
coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2];
n_coord_ = 1;
this->r() += TINY_BIT * this->u();
// Couldn't find next cell anywhere! This probably means there is an actual
// undefined region in the geometry.
if (!find_cell(this, false)) {
this->mark_as_lost("After particle " + std::to_string(id) +
this->mark_as_lost("After particle " + std::to_string(id_) +
" crossed surface " + std::to_string(surf->id_) +
" it could not be located in any cell and it did not leak.");
return;
@ -631,7 +588,7 @@ Particle::mark_as_lost(const char* message)
write_restart();
// Increment number of lost particles
alive = false;
alive_ = false;
#pragma omp atomic
simulation::n_lost_particles += 1;
@ -655,7 +612,7 @@ Particle::write_restart() const
// Set up file name
std::stringstream filename;
filename << settings::path_output << "particle_" << simulation::current_batch
<< '_' << id << ".h5";
<< '_' << id_ << ".h5";
#pragma omp critical (WriteParticleRestart)
{
@ -686,15 +643,14 @@ Particle::write_restart() const
write_dataset(file_id, "run_mode", "particle restart");
break;
}
write_dataset(file_id, "id", id);
write_dataset(file_id, "type", type);
write_dataset(file_id, "id", id_);
write_dataset(file_id, "type", static_cast<int>(type_));
int64_t i = simulation::current_work;
write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
hsize_t dims[] {3};
write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false);
write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false);
write_dataset(file_id, "xyz", simulation::source_bank[i-1].r);
write_dataset(file_id, "uvw", simulation::source_bank[i-1].u);
// Close file
file_close(file_id);

38
src/particle_restart.cpp
View file

@ -39,29 +39,28 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
} else if (mode == "fixed source") {
previous_run_mode = RUN_MODE_FIXEDSOURCE;
}
read_dataset(file_id, "id", p.id);
read_dataset(file_id, "type", p.type);
read_dataset(file_id, "weight", p.wgt);
read_dataset(file_id, "energy", p.E);
std::array<double, 3> x;
read_dataset(file_id, "xyz", x);
std::copy(x.data(), x.data() + 3, p.coord[0].xyz);
read_dataset(file_id, "uvw", x);
std::copy(x.data(), x.data() + 3, p.coord[0].uvw);
read_dataset(file_id, "id", p.id_);
int type;
read_dataset(file_id, "type", type);
p.type_ = static_cast<Particle::Type>(type);
read_dataset(file_id, "weight", p.wgt_);
read_dataset(file_id, "energy", p.E_);
read_dataset(file_id, "xyz", p.r());
read_dataset(file_id, "uvw", p.u());
// Set energy group and average energy in multi-group mode
if (!settings::run_CE) {
p.g = p.E;
p.E = data::energy_bin_avg[p.g - 1];
p.g_ = p.E_;
p.E_ = data::energy_bin_avg[p.g_ - 1];
}
// Set particle last attributes
p.last_wgt = p.wgt;
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current);
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz);
std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw);
p.last_E = p.E;
p.last_g = p.g;
p.wgt_last_ = p.wgt_;
p.r_last_current_ = p.r();
p.r_last_ = p.r();
p.u_last_ = p.u();
p.E_last_ = p.E_;
p.g_last_ = p.g_;
// Close hdf5 file
file_close(file_id);
@ -81,7 +80,6 @@ void run_particle_restart()
// Initialize the particle to be tracked
Particle p;
p.initialize();
// Read in the restart information
int previous_run_mode;
@ -94,10 +92,10 @@ void run_particle_restart()
int64_t particle_seed;
switch (previous_run_mode) {
case RUN_MODE_EIGENVALUE:
particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id;
particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id_;
break;
case RUN_MODE_FIXEDSOURCE:
particle_seed = p.id;
particle_seed = p.id_;
break;
}
set_particle_seed(particle_seed);

27
src/photon.cpp
View file

@ -224,7 +224,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
}
// Truncate the bremsstrahlung data at the cutoff energy
int photon = static_cast<int>(ParticleType::photon);
int photon = static_cast<int>(Particle::Type::photon);
const auto& E {electron_energy};
double cutoff = settings::energy_cutoff[photon];
if (cutoff > E(0)) {
@ -653,13 +653,12 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
if (shell.n_transitions == 0) {
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<double, 3> uvw;
uvw[0] = mu;
uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi);
uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
double E = shell.binding_energy;
int photon = static_cast<int>(ParticleType::photon);
p.create_secondary(uvw.data(), E, photon, true);
p.create_secondary(u, E, Particle::Type::photon);
return;
}
@ -679,10 +678,10 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Sample angle isotropically
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<double, 3> uvw;
uvw[0] = mu;
uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi);
uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Get the transition energy
double E = shell.transition_energy(i_transition);
@ -691,8 +690,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Non-radiative transition -- Auger/Coster-Kronig effect
// Create auger electron
int electron = static_cast<int>(ParticleType::electron);
p.create_secondary(uvw.data(), E, electron, true);
p.create_secondary(u, E, Particle::Type::electron);
// Fill hole left by emitted auger electron
int i_hole = shell_map_.at(secondary);
@ -702,8 +700,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Radiative transition -- get X-ray energy
// Create fluorescent photon
int photon = static_cast<int>(ParticleType::photon);
p.create_secondary(uvw.data(), E, photon, true);
p.create_secondary(u, E, Particle::Type::photon);
}
// Fill hole created by electron transitioning to the photoelectron hole

296
src/physics.cpp
View file

@ -33,40 +33,41 @@ namespace openmc {
void collision(Particle* p)
{
// Add to collision counter for particle
++(p->n_collision);
++(p->n_collision_);
// Sample reaction for the material the particle is in
switch (static_cast<ParticleType>(p->type)) {
case ParticleType::neutron:
switch (static_cast<Particle::Type>(p->type_)) {
case Particle::Type::neutron:
sample_neutron_reaction(p);
break;
case ParticleType::photon:
case Particle::Type::photon:
sample_photon_reaction(p);
break;
case ParticleType::electron:
case Particle::Type::electron:
sample_electron_reaction(p);
break;
case ParticleType::positron:
case Particle::Type::positron:
sample_positron_reaction(p);
break;
}
// Kill particle if energy falls below cutoff
if (p->E < settings::energy_cutoff[p->type]) {
p->alive = false;
p->wgt = 0.0;
p->last_wgt = 0.0;
int type = static_cast<int>(p->type_);
if (p->E_ < settings::energy_cutoff[type]) {
p->alive_ = false;
p->wgt_ = 0.0;
p->wgt_last_ = 0.0;
}
// Display information about collision
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
if (static_cast<ParticleType>(p->type) == ParticleType::neutron) {
msg << " " << reaction_name(p->event_MT) << " with " <<
data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV.";
if (p->type_ == Particle::Type::neutron) {
msg << " " << reaction_name(p->event_mt_) << " with " <<
data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV.";
} else {
msg << " " << reaction_name(p->event_MT) << " with " <<
data::elements[p->event_nuclide].name_ << ". Energy = " << p->E << " eV.";
msg << " " << reaction_name(p->event_mt_) << " with " <<
data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV.";
}
write_message(msg, 1);
}
@ -78,7 +79,7 @@ void sample_neutron_reaction(Particle* p)
int i_nuclide = sample_nuclide(p);
// Save which nuclide particle had collision with
p->event_nuclide = i_nuclide;
p->event_nuclide_ = i_nuclide;
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
@ -88,14 +89,14 @@ void sample_neutron_reaction(Particle* p)
const auto& nuc {data::nuclides[i_nuclide]};
if (nuc->fissionable_) {
Reaction* rx = sample_fission(i_nuclide, p->E);
Reaction* rx = sample_fission(i_nuclide, p->E_);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(),
&simulation::n_bank, simulation::fission_bank.size());
} else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
settings::create_fission_neutrons) {
create_fission_sites(p, i_nuclide, rx, p->secondary_bank,
&p->n_secondary, MAX_SECONDARY);
create_fission_sites(p, i_nuclide, rx, p->secondary_bank_,
&p->n_secondary_, MAX_SECONDARY);
}
}
@ -112,16 +113,16 @@ void sample_neutron_reaction(Particle* p)
if (simulation::micro_xs[i_nuclide].absorption > 0.0) {
absorption(p, i_nuclide);
} else {
p->absorb_wgt = 0.0;
p->wgt_absorb_ = 0.0;
}
if (!p->alive) return;
if (!p->alive_) return;
// Sample a scattering reaction and determine the secondary energy of the
// exiting neutron
scatter(p, i_nuclide);
// Advance URR seed stream 'N' times after energy changes
if (p->E != p->last_E) {
if (p->E_ != p->E_last_) {
prn_set_stream(STREAM_URR_PTABLE);
advance_prn_seed(data::nuclides.size());
prn_set_stream(STREAM_TRACKING);
@ -130,13 +131,13 @@ void sample_neutron_reaction(Particle* p)
// Play russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p->alive) return;
if (!p->alive_) return;
}
}
void
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array,
int64_t* size_bank, int64_t bank_capacity)
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity)
{
// TODO: Heat generation from fission
@ -145,7 +146,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt / simulation::keff * weight * simulation::micro_xs[
double nu_t = p->wgt_ / simulation::keff * weight * simulation::micro_xs[
i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].total;
// Sample the number of neutrons produced
@ -177,28 +178,26 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission = true;
p->fission_ = true;
for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) {
// Bank source neutrons by copying the particle data
bank_array[i].xyz[0] = p->coord[0].xyz[0];
bank_array[i].xyz[1] = p->coord[0].xyz[1];
bank_array[i].xyz[2] = p->coord[0].xyz[2];
bank_array[i].r = p->r();
// Set that the bank particle is a neutron
bank_array[i].particle = static_cast<int>(ParticleType::neutron);
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E, &bank_array[i]);
sample_fission_neutron(i_nuclide, rx, p->E_, &bank_array[i]);
// Set the delayed group on the particle as well
p->delayed_group = bank_array[i].delayed_group;
p->delayed_group_ = bank_array[i].delayed_group;
// Increment the number of neutrons born delayed
if (p->delayed_group > 0) {
nu_d[p->delayed_group-1]++;
if (p->delayed_group_ > 0) {
nu_d[p->delayed_group_-1]++;
}
}
@ -206,10 +205,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_
*size_bank = std::min(*size_bank + nu, bank_capacity);
// Store the total weight banked for analog fission tallies
p->n_bank = nu;
p->wgt_bank = nu / weight;
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
p->n_delayed_bank[d] = nu_d[d];
p->n_delayed_bank_[d] = nu_d[d];
}
}
@ -218,21 +217,21 @@ void sample_photon_reaction(Particle* p)
// Kill photon if below energy cutoff -- an extra check is made here because
// photons with energy below the cutoff may have been produced by neutrons
// reactions or atomic relaxation
int photon = static_cast<int>(ParticleType::photon);
if (p->E < settings::energy_cutoff[photon]) {
p->E = 0.0;
p->alive = false;
int photon = static_cast<int>(Particle::Type::photon);
if (p->E_ < settings::energy_cutoff[photon]) {
p->E_ = 0.0;
p->alive_ = false;
return;
}
// Sample element within material
int i_element = sample_element(p);
p->event_nuclide = i_element;
p->event_nuclide_ = i_element;
const auto& micro {simulation::micro_photon_xs[i_element]};
const auto& element {data::elements[i_element]};
// Calculate photon energy over electron rest mass equivalent
double alpha = p->E/MASS_ELECTRON_EV;
double alpha = p->E_/MASS_ELECTRON_EV;
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
@ -243,8 +242,8 @@ void sample_photon_reaction(Particle* p)
prob += micro.coherent;
if (prob > cutoff) {
double mu = element.rayleigh_scatter(alpha);
rotate_angle_c(p->coord[0].uvw, mu, nullptr);
p->event_MT = COHERENT;
p->u() = rotate_angle(p->u(), mu, nullptr);
p->event_mt_ = COHERENT;
return;
}
@ -269,11 +268,8 @@ void sample_photon_reaction(Particle* p)
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
double phi = 2.0*PI*prn();
double uvw[3];
std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_electron, &phi);
int electron = static_cast<int>(ParticleType::electron);
p->create_secondary(uvw, E_electron, electron, true);
Direction u = rotate_angle(p->u(), mu_electron, &phi);
p->create_secondary(u, E_electron, Particle::Type::electron);
// TODO: Compton subshell data does not match atomic relaxation data
// Allow electrons to fill orbital and produce auger electrons
@ -284,9 +280,9 @@ void sample_photon_reaction(Particle* p)
}
phi += PI;
p->E = alpha_out*MASS_ELECTRON_EV;
rotate_angle_c(p->coord[0].uvw, mu, &phi);
p->event_MT = INCOHERENT;
p->E_ = alpha_out*MASS_ELECTRON_EV;
p->u() = rotate_angle(p->u(), mu, &phi);
p->event_mt_ = INCOHERENT;
return;
}
@ -309,7 +305,7 @@ void sample_photon_reaction(Particle* p)
prob += xs;
if (prob > cutoff) {
double E_electron = p->E - shell.binding_energy;
double E_electron = p->E_ - shell.binding_energy;
// Sample mu using non-relativistic Sauter distribution.
// See Eqns 3.19 and 3.20 in "Implementing a photon physics
@ -326,21 +322,20 @@ void sample_photon_reaction(Particle* p)
}
double phi = 2.0*PI*prn();
std::array<double, 3> uvw;
uvw[0] = mu;
uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi);
uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create secondary electron
int electron = static_cast<int>(ParticleType::electron);
p->create_secondary(uvw.data(), E_electron, electron, true);
p->create_secondary(u, E_electron, Particle::Type::electron);
// Allow electrons to fill orbital and produce auger electrons
// and fluorescent photons
element.atomic_relaxation(shell, *p);
p->event_MT = 533 + shell.index_subshell;
p->alive = false;
p->E = 0.0;
p->event_mt_ = 533 + shell.index_subshell;
p->alive_ = false;
p->E_ = 0.0;
return;
}
}
@ -356,21 +351,16 @@ void sample_photon_reaction(Particle* p)
&mu_electron, &mu_positron);
// Create secondary electron
double uvw[3];
std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_electron, nullptr);
int electron = static_cast<int>(ParticleType::electron);
p->create_secondary(uvw, E_electron, electron, true);
Direction u = rotate_angle(p->u(), mu_electron, nullptr);
p->create_secondary(u, E_electron, Particle::Type::electron);
// Create secondary positron
std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_positron, nullptr);
int positron = static_cast<int>(ParticleType::positron);
p->create_secondary(uvw, E_positron, positron, true);
u = rotate_angle(p->u(), mu_positron, nullptr);
p->create_secondary(u, E_positron, Particle::Type::positron);
p->event_MT = PAIR_PROD;
p->alive = false;
p->E = 0.0;
p->event_mt_ = PAIR_PROD;
p->alive_ = false;
p->E_ = 0.0;
}
}
@ -383,8 +373,8 @@ void sample_electron_reaction(Particle* p)
thick_target_bremsstrahlung(*p, &E_lost);
}
p->E = 0.0;
p->alive = false;
p->E_ = 0.0;
p->alive_ = false;
}
void sample_positron_reaction(Particle* p)
@ -399,22 +389,17 @@ void sample_positron_reaction(Particle* p)
// Sample angle isotropically
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<double, 3> uvw;
uvw[0] = mu;
uvw[1] = std::sqrt(1.0 - mu*mu)*std::cos(phi);
uvw[2] = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create annihilation photon pair traveling in opposite directions
int photon = static_cast<int>(ParticleType::photon);
p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true);
p->create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon);
p->create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon);
uvw[0] = -uvw[0];
uvw[1] = -uvw[1];
uvw[2] = -uvw[2];
p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true);
p->E = 0.0;
p->alive = false;
p->E_ = 0.0;
p->alive_ = false;
}
int sample_nuclide(const Particle* p)
@ -423,7 +408,7 @@ int sample_nuclide(const Particle* p)
double cutoff = prn() * simulation::material_xs.total;
// Get pointers to nuclide/density arrays
const auto& mat {model::materials[p->material]};
const auto& mat {model::materials[p->material_]};
int n = mat->nuclide_.size();
double prob = 0.0;
@ -448,7 +433,7 @@ int sample_element(Particle* p)
double cutoff = prn() * simulation::material_xs.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material]};
const auto& mat {model::materials[p->material_]};
int n = mat->nuclide_.size();
int i = 0;
@ -541,7 +526,7 @@ void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product)
+ f*(rx->xs_[i_temp].value[i_grid - threshold + 1]));
for (int j = 0; j < rx->products_.size(); ++j) {
if (rx->products_[j].particle_ == ParticleType::photon) {
if (rx->products_[j].particle_ == Particle::Type::photon) {
// add to cumulative probability
prob += (*rx->products_[j].yield_)(E) * xs;
@ -557,16 +542,16 @@ void absorption(Particle* p, int i_nuclide)
{
if (settings::survival_biasing) {
// Determine weight absorbed in survival biasing
p->absorb_wgt = p->wgt * simulation::micro_xs[i_nuclide].absorption /
p->wgt_absorb_ = p->wgt_ * simulation::micro_xs[i_nuclide].absorption /
simulation::micro_xs[i_nuclide].total;
// Adjust weight of particle by probability of absorption
p->wgt -= p->absorb_wgt;
p->last_wgt = p->wgt;
p->wgt_ -= p->wgt_absorb_;
p->wgt_last_ = p->wgt_;
// Score implicit absorption estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
global_tally_absorption += p->absorb_wgt * simulation::micro_xs[
global_tally_absorption += p->wgt_absorb_ * simulation::micro_xs[
i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption;
}
} else {
@ -575,13 +560,13 @@ void absorption(Particle* p, int i_nuclide)
prn() * simulation::micro_xs[i_nuclide].total) {
// Score absorption estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
global_tally_absorption += p->wgt * simulation::micro_xs[
global_tally_absorption += p->wgt_ * simulation::micro_xs[
i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption;
}
p->alive = false;
p->event = EVENT_ABSORB;
p->event_MT = N_DISAPPEAR;
p->alive_ = false;
p->event_ = EVENT_ABSORB;
p->event_mt_ = N_DISAPPEAR;
}
}
}
@ -589,7 +574,7 @@ void absorption(Particle* p, int i_nuclide)
void scatter(Particle* p, int i_nuclide)
{
// copy incoming direction
Direction u_old {p->coord[0].uvw};
Direction u_old {p->u()};
// Get pointer to nuclide and grid index/interpolation factor
const auto& nuc {data::nuclides[i_nuclide]};
@ -614,13 +599,13 @@ void scatter(Particle* p, int i_nuclide)
// NON-S(A,B) ELASTIC SCATTERING
// Determine temperature
double kT = nuc->multipole_ ? p->sqrtkT*p->sqrtkT : nuc->kTs_[i_temp];
double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT,
&p->E, p->coord[0].uvw, &p->mu, &p->wgt);
elastic_scatter(i_nuclide, *nuc->reactions_[0], kT,
p->E_, p->u(), p->mu_);
p->event_MT = ELASTIC;
p->event_mt_ = ELASTIC;
sampled = true;
}
@ -629,9 +614,9 @@ void scatter(Particle* p, int i_nuclide)
// =======================================================================
// S(A,B) SCATTERING
sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu);
sab_scatter(i_nuclide, micro.index_sab, p->E_, p->u(), p->mu_);
p->event_MT = ELASTIC;
p->event_mt_ = ELASTIC;
sampled = true;
}
@ -663,53 +648,47 @@ void scatter(Particle* p, int i_nuclide)
// Perform collision physics for inelastic scattering
const auto& rx {nuc->reactions_[i]};
inelastic_scatter(nuc.get(), rx.get(), p);
p->event_MT = rx->mt_;
p->event_mt_ = rx->mt_;
}
// Set event component
p->event = EVENT_SCATTER;
p->event_ = EVENT_SCATTER;
// Sample new outgoing angle for isotropic-in-lab scattering
const auto& mat {model::materials[p->material]};
const auto& mat {model::materials[p->material_]};
if (!mat->p0_.empty()) {
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
if (mat->p0_[i_nuc_mat]) {
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
Direction u_new;
u_new.x = mu;
u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p->mu = u_old.dot(u_new);
// Change direction of particle
p->coord[0].uvw[0] = u_new.x;
p->coord[0].uvw[1] = u_new.y;
p->coord[0].uvw[2] = u_new.z;
p->u().x = mu;
p->u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
p->u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p->mu_ = u_old.dot(p->u());
}
}
}
void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
double* uvw, double* mu_lab, double* wgt)
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E,
Direction& u, double& mu_lab)
{
// get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide]};
double vel = std::sqrt(*E);
double vel = std::sqrt(E);
double awr = nuc->awr_;
// Neutron velocity in LAB
Direction u {uvw};
Direction v_n = vel*u;
// Sample velocity of target nucleus
Direction v_t {};
if (!simulation::micro_xs[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), *E, u, v_n,
simulation::micro_xs[i_nuclide].elastic, kT, wgt);
v_t = sample_target_velocity(nuc.get(), E, u, v_n,
simulation::micro_xs[i_nuclide].elastic, kT);
}
// Velocity of center-of-mass
@ -724,10 +703,10 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
// Sample scattering angle, checking if it is an ncorrelated angle-energy
// distribution
double mu_cm;
auto& d = rx->products_[0].distribution_[0];
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) {
mu_cm = d_->angle().sample(*E);
mu_cm = d_->angle().sample(E);
} else {
mu_cm = 2.0*prn() - 1.0;
}
@ -743,26 +722,23 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
// Transform back to LAB frame
v_n += v_cm;
*E = v_n.dot(v_n);
vel = std::sqrt(*E);
E = v_n.dot(v_n);
vel = std::sqrt(E);
// compute cosine of scattering angle in LAB frame by taking dot product of
// neutron's pre- and post-collision angle
*mu_lab = u.dot(v_n) / vel;
mu_lab = u.dot(v_n) / vel;
// Set energy and direction of particle in LAB frame
u = v_n / vel;
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
// Because of floating-point roundoff, it may be possible for mu_lab to be
// outside of the range [-1,1). In these cases, we just set mu_lab to exactly
// -1 or 1
if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab);
if (std::abs(mu_lab) > 1.0) mu_lab = std::copysign(1.0, mu_lab);
}
void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu)
void sab_scatter(int i_nuclide, int i_sab, double& E, Direction& u, double& mu)
{
// Determine temperature index
const auto& micro {simulation::micro_xs[i_nuclide]};
@ -770,15 +746,15 @@ void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu)
// Sample energy and angle
double E_out;
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, *E, &E_out, mu);
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, E, &E_out, &mu);
// Set energy to outgoing, change direction of particle
*E = E_out;
rotate_angle_c(uvw, *mu, nullptr);
E = E_out;
u = rotate_angle(u, mu, nullptr);
}
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_neut, double xs_eff, double kT, double* wgt)
Direction v_neut, double xs_eff, double kT)
{
// check if nuclide is a resonant scatterer
ResScatMethod sampling_method;
@ -977,7 +953,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT)
return vt * rotate_angle(u, mu, nullptr);
}
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site)
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site)
{
// Sample cosine of angle -- fission neutrons are always emitted
// isotropically. Sometimes in ACE data, fission reactions actually have
@ -987,9 +963,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
// Sample azimuthal angle uniformly in [0,2*pi)
double phi = 2.0*PI*prn();
site->uvw[0] = mu;
site->uvw[1] = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->uvw[2] = std::sqrt(1.0 - mu*mu) * std::sin(phi);
site->u.x = mu;
site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
// Determine total nu, delayed nu, and delayed neutron fraction
const auto& nuc {data::nuclides[i_nuclide]};
@ -1029,7 +1005,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
rx->products_[group].sample(E_in, site->E, mu);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(ParticleType::neutron);
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
if (site->E < data::energy_max[neutron]) break;
// check for large number of resamples
@ -1054,7 +1030,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
rx->products_[0].sample(E_in, site->E, mu);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(ParticleType::neutron);
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
if (site->E < data::energy_max[neutron]) break;
// check for large number of resamples
@ -1071,7 +1047,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
{
// copy energy of neutron
double E_in = p->E;
double E_in = p->E_;
// sample outgoing energy and scattering cosine
double E;
@ -1098,30 +1074,29 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
// Set outgoing energy and scattering angle
p->E = E;
p->mu = mu;
p->E_ = E;
p->mu_ = mu;
// change direction of particle
rotate_angle_c(p->coord[0].uvw, mu, nullptr);
p->u() = rotate_angle(p->u(), mu, nullptr);
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
if (std::floor(yield) == yield) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
int neutron = static_cast<int>(ParticleType::neutron);
p->create_secondary(p->coord[0].uvw, p->E, neutron, true);
p->create_secondary(p->u(), p->E_, Particle::Type::neutron);
}
} else {
// Otherwise, change weight of particle based on yield
p->wgt *= yield;
p->wgt_ *= yield;
}
}
void sample_secondary_photons(Particle* p, int i_nuclide)
{
// Sample the number of photons produced
double y_t = p->wgt * simulation::micro_xs[i_nuclide].photon_prod /
double y_t = p->wgt_ * simulation::micro_xs[i_nuclide].photon_prod /
simulation::micro_xs[i_nuclide].total;
int y = static_cast<int>(y_t);
if (prn() <= y_t - y) ++y;
@ -1131,22 +1106,19 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
// Sample the reaction and product
int i_rx;
int i_product;
sample_photon_product(i_nuclide, p->E, &i_rx, &i_product);
sample_photon_product(i_nuclide, p->E_, &i_rx, &i_product);
// Sample the outgoing energy and angle
auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
double E;
double mu;
rx->products_[i_product].sample(p->E, E, mu);
rx->products_[i_product].sample(p->E_, E, mu);
// Sample the new direction
double uvw[3];
std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu, nullptr);
Direction u = rotate_angle(p->u(), mu, nullptr);
// Create the secondary photon
int photon = static_cast<int>(ParticleType::photon);
p->create_secondary(uvw, E, photon, true);
p->create_secondary(u, E, Particle::Type::photon);
}
}

14
src/physics_common.cpp
View file

@ -11,14 +11,14 @@ namespace openmc {
void russian_roulette(Particle* p)
{
if (p->wgt < settings::weight_cutoff) {
if (prn() < p->wgt / settings::weight_survive) {
p->wgt = settings::weight_survive;
p->last_wgt = p->wgt;
if (p->wgt_ < settings::weight_cutoff) {
if (prn() < p->wgt_ / settings::weight_survive) {
p->wgt_ = settings::weight_survive;
p->wgt_last_ = p->wgt_;
} else {
p->wgt = 0.;
p->last_wgt = 0.;
p->alive = false;
p->wgt_ = 0.;
p->wgt_last_ = 0.;
p->alive_ = false;
}
}
}

76
src/physics_mg.cpp
View file

@ -25,7 +25,7 @@ void
collision_mg(Particle* p)
{
// Add to the collision counter for the particle
p->n_collision++;
p->n_collision_++;
// Sample the reaction type
sample_reaction(p);
@ -33,7 +33,7 @@ collision_mg(Particle* p)
// Display information about collision
if ((settings::verbosity >= 10) || (simulation::trace)) {
std::stringstream msg;
msg << " Energy Group = " << p->g;
msg << " Energy Group = " << p->g_;
write_message(msg, 1);
}
}
@ -46,14 +46,14 @@ sample_reaction(Particle* p)
// change when sampling fission sites. The following block handles all
// absorption (including fission)
if (model::materials[p->material]->fissionable_) {
if (model::materials[p->material_]->fissionable_) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
create_fission_sites(
p, simulation::fission_bank.data(), &simulation::n_bank,
simulation::fission_bank.size());
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_),
MAX_SECONDARY);
}
}
@ -63,9 +63,9 @@ sample_reaction(Particle* p)
if (simulation::material_xs.absorption > 0.) {
absorption(p);
} else {
p->absorb_wgt = 0.;
p->wgt_absorb_ = 0.;
}
if (!p->alive) return;
if (!p->alive_) return;
// Sample a scattering event to determine the energy of the exiting neutron
scatter(p);
@ -73,7 +73,7 @@ sample_reaction(Particle* p)
// Play Russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p->alive) return;
if (!p->alive_) return;
}
}
@ -82,27 +82,27 @@ scatter(Particle* p)
{
// Adjust indices for Fortran to C++ indexing
// TODO: Remove when no longer needed
int gin = p->last_g - 1;
int gout = p->g - 1;
int i_mat = p->material;
data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
int gin = p->g_last_ - 1;
int gout = p->g_ - 1;
int i_mat = p->material_;
data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_);
// Adjust return value for fortran indexing
// TODO: Remove when no longer needed
p->g = gout + 1;
p->g_ = gout + 1;
// Rotate the angle
rotate_angle_c(p->coord[0].uvw, p->mu, nullptr);
p->u() = rotate_angle(p->u(), p->mu_, nullptr);
// Update energy value for downstream compatability (in tallying)
p->E = data::energy_bin_avg[gout];
p->E_ = data::energy_bin_avg[gout];
// Set event component
p->event = EVENT_SCATTER;
p->event_ = EVENT_SCATTER;
}
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size)
{
// TODO: Heat generation from fission
@ -112,7 +112,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt / simulation::keff * weight *
double nu_t = p->wgt_ / simulation::keff * weight *
simulation::material_xs.nu_fission / simulation::material_xs.total;
// Sample the number of neutrons produced
@ -148,16 +148,14 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
// group.
double nu_d[MAX_DELAYED_GROUPS] = {0.};
p->fission = true;
p->fission_ = true;
for (size_t i = static_cast<size_t>(*size_bank);
i < static_cast<size_t>(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<int>(ParticleType::neutron);
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
@ -168,23 +166,23 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn();
bank_array[i].uvw[0] = mu;
bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi);
bank_array[i].u.x = mu;
bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
// Sample secondary energy distribution for the fission reaction and set
// the energy in the fission bank
int dg;
int gout;
data::macro_xs[p->material].sample_fission_energy(p->g - 1, dg, gout);
bank_array[i].E = static_cast<double>(gout + 1);
data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
bank_array[i].E = gout + 1;
bank_array[i].delayed_group = dg + 1;
// Set the delayed group on the particle as well
p->delayed_group = dg + 1;
p->delayed_group_ = dg + 1;
// Increment the number of neutrons born delayed
if (p->delayed_group > 0) {
if (p->delayed_group_ > 0) {
nu_d[dg]++;
}
}
@ -193,10 +191,10 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
*size_bank = std::min(*size_bank + nu, bank_array_size);
// Store the total weight banked for analog fission tallies
p->n_bank = nu;
p->wgt_bank = nu / weight;
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
p->n_delayed_bank[d] = nu_d[d];
p->n_delayed_bank_[d] = nu_d[d];
}
}
@ -205,26 +203,26 @@ absorption(Particle* p)
{
if (settings::survival_biasing) {
// Determine weight absorbed in survival biasing
p->absorb_wgt = p->wgt *
p->wgt_absorb_ = p->wgt_ *
simulation::material_xs.absorption / simulation::material_xs.total;
// Adjust weight of particle by the probability of absorption
p->wgt -= p->absorb_wgt;
p->last_wgt = p->wgt;
p->wgt_ -= p->wgt_absorb_;
p->wgt_last_ = p->wgt_;
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->absorb_wgt *
global_tally_absorption += p->wgt_absorb_ *
simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
} else {
if (simulation::material_xs.absorption >
prn() * simulation::material_xs.total) {
#pragma omp atomic
global_tally_absorption += p->wgt * simulation::material_xs.nu_fission /
global_tally_absorption += p->wgt_ * simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
p->alive = false;
p->event = EVENT_ABSORB;
p->alive_ = false;
p->event_ = EVENT_ABSORB;
}
}

140
src/plot.cpp
View file

@ -106,49 +106,48 @@ void create_ppm(Plot pl)
data.resize({width, height});
int in_i, out_i;
double xyz[3];
Position r;
switch(pl.basis_) {
case PlotBasis::xy :
in_i = 0;
out_i = 1;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1] + pl.width_[1] / 2.;
xyz[2] = pl.origin_[2];
r.x = pl.origin_[0] - pl.width_[0] / 2.;
r.y = pl.origin_[1] + pl.width_[1] / 2.;
r.z = pl.origin_[2];
break;
case PlotBasis::xz :
in_i = 0;
out_i = 2;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1];
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
r.x = pl.origin_[0] - pl.width_[0] / 2.;
r.y = pl.origin_[1];
r.z = pl.origin_[2] + pl.width_[1] / 2.;
break;
case PlotBasis::yz :
in_i = 1;
out_i = 2;
xyz[0] = pl.origin_[0];
xyz[1] = pl.origin_[1] - pl.width_[0] / 2.;
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
r.x = pl.origin_[0];
r.y = pl.origin_[1] - pl.width_[0] / 2.;
r.z = pl.origin_[2] + pl.width_[1] / 2.;
break;
}
double dir[3] = {0.5, 0.5, 0.5};
Direction u {0.5, 0.5, 0.5};
#pragma omp parallel
{
Particle p;
p.initialize();
std::copy(xyz, xyz+3, p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = model::root_universe;
p.r() = r;
p.u() = u;
p.coord_[0].universe = model::root_universe;
#pragma omp for
for (int y = 0; y < height; y++) {
p.coord[0].xyz[out_i] = xyz[out_i] - out_pixel * y;
p.r()[out_i] = r[out_i] - out_pixel * y;
for (int x = 0; x < width; x++) {
// local variables
RGBColor rgb;
int id;
p.coord[0].xyz[in_i] = xyz[in_i] + in_pixel * x;
p.r()[in_i] = r[in_i] + in_pixel * x;
position_rgb(p, pl, rgb, id);
data(x,y) = rgb;
}
@ -309,11 +308,9 @@ void
Plot::set_origin(pugi::xml_node plot_node)
{
// Copy plotting origin
std::vector<double> pl_origin = get_node_array<double>(plot_node, "origin");
auto pl_origin = get_node_array<double>(plot_node, "origin");
if (pl_origin.size() == 3) {
origin_[0] = pl_origin[0];
origin_[1] = pl_origin[1];
origin_[2] = pl_origin[2];
origin_ = pl_origin;
} else {
std::stringstream err_msg;
err_msg << "Origin must be length 3 in plot "
@ -329,8 +326,8 @@ Plot::set_width(pugi::xml_node plot_node)
std::vector<double> pl_width = get_node_array<double>(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 << "<width> must be length 2 in slice plot "
@ -340,9 +337,7 @@ Plot::set_width(pugi::xml_node plot_node)
} else if (PlotType::voxel == type_) {
if (pl_width.size() == 3) {
pl_width = get_node_array<double>(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 << "<width> must be length 3 in voxel plot "
@ -648,11 +643,11 @@ index_meshlines_mesh_(-1)
void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
{
p.n_coord = 1;
p.n_coord_ = 1;
bool found_cell = find_cell(&p, 0);
int j = p.n_coord - 1;
int j = p.n_coord_ - 1;
if (settings::check_overlaps) {check_cell_overlap(&p);}
@ -666,23 +661,23 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id)
} else {
if (PlotColorBy::mats == pl.color_by_) {
// Assign color based on material
Cell* c = model::cells[p.coord[j].cell];
const auto& c = model::cells[p.coord_[j].cell];
if (c->type_ == FILL_UNIVERSE) {
// If we stopped on a middle universe level, treat as if not found
rgb = pl.not_found_;
id = -1;
} else if (p.material == MATERIAL_VOID) {
} else if (p.material_ == MATERIAL_VOID) {
// By default, color void cells white
rgb = WHITE;
id = -1;
} else {
rgb = pl.colors_[p.material];
id = model::materials[p.material]->id_;
rgb = pl.colors_[p.material_];
id = model::materials[p.material_]->id_;
}
} else if (PlotColorBy::cells == pl.color_by_) {
// Assign color based on cell
rgb = pl.colors_[p.coord[j].cell];
id = model::cells[p.coord[j].cell]->id_;
rgb = pl.colors_[p.coord_[j].cell];
id = model::cells[p.coord_[j].cell]->id_;
}
} // endif found_cell
}
@ -746,36 +741,27 @@ void draw_mesh_lines(Plot pl, ImageData& data)
break;
}
double xyz_ll_plot[3], xyz_ur_plot[3];
xyz_ll_plot[0] = pl.origin_[0];
xyz_ll_plot[1] = pl.origin_[1];
xyz_ll_plot[2] = pl.origin_[2];
Position ll_plot {pl.origin_};
Position ur_plot {pl.origin_};
xyz_ur_plot[0] = pl.origin_[0];
xyz_ur_plot[1] = pl.origin_[1];
xyz_ur_plot[2] = pl.origin_[2];
ll_plot[outer] -= pl.width_[0] / 2.;
ll_plot[inner] -= pl.width_[1] / 2.;
ur_plot[outer] += pl.width_[0] / 2.;
ur_plot[inner] += pl.width_[1] / 2.;
xyz_ll_plot[outer] = pl.origin_[outer] - pl.width_[0] / 2.;
xyz_ll_plot[inner] = pl.origin_[inner] - pl.width_[1] / 2.;
xyz_ur_plot[outer] = pl.origin_[outer] + pl.width_[0] / 2.;
xyz_ur_plot[inner] = pl.origin_[inner] + pl.width_[1] / 2.;
int width[3];
width[0] = xyz_ur_plot[0] - xyz_ll_plot[0];
width[1] = xyz_ur_plot[1] - xyz_ll_plot[1];
width[2] = xyz_ur_plot[2] - xyz_ll_plot[2];
Position width = ur_plot - ll_plot;
auto& m = model::meshes[pl.index_meshlines_mesh_];
int ijk_ll[3], ijk_ur[3];
bool in_mesh;
m->get_indices(Position(xyz_ll_plot), &(ijk_ll[0]), &in_mesh);
m->get_indices(Position(xyz_ur_plot), &(ijk_ur[0]), &in_mesh);
m->get_indices(ll_plot, &(ijk_ll[0]), &in_mesh);
m->get_indices(ur_plot, &(ijk_ur[0]), &in_mesh);
// Fortran/C++ index correction
ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++;
double xyz_ll[3], xyz_ur[3];
Position r_ll, r_ur;
// sweep through all meshbins on this plane and draw borders
for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) {
for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) {
@ -783,20 +769,20 @@ void draw_mesh_lines(Plot pl, ImageData& data)
if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) {
int outrange[3], inrange[3];
// get xyz's of lower left and upper right of this mesh cell
xyz_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
xyz_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
xyz_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
xyz_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
r_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
r_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
r_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
r_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
// map the xyz ranges to pixel ranges
double frac = (xyz_ll[outer] - xyz_ll_plot[outer]) / width[outer];
double frac = (r_ll[outer] - ll_plot[outer]) / width[outer];
outrange[0] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[outer] - xyz_ll_plot[outer]) / width[outer];
frac = (r_ur[outer] - ll_plot[outer]) / width[outer];
outrange[1] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[inner] - xyz_ll_plot[inner]) / width[inner];
frac = (r_ur[inner] - ll_plot[inner]) / width[inner];
inrange[0] = int((1. - frac) * (double)pl.pixels_[1]);
frac = (xyz_ll[inner] - xyz_ll_plot[inner]) / width[inner];
frac = (r_ll[inner] - ll_plot[inner]) / width[inner];
inrange[1] = int((1. - frac) * (double)pl.pixels_[1]);
// draw lines
@ -846,18 +832,14 @@ void create_voxel(Plot pl)
vox[2] = pl.width_[2]/(double)pl.pixels_[2];
// initial particle position
std::array<double, 3> ll;
ll[0] = pl.origin_[0] - pl.width_[0] / 2.;
ll[1] = pl.origin_[1] - pl.width_[1] / 2.;
ll[2] = pl.origin_[2] - pl.width_[2] / 2.;
Position ll = pl.origin_ - pl.width_ / 2.;
// allocate and initialize particle
double dir[3] = {0.5, 0.5, 0.5};
Direction u {0.5, 0.5, 0.5};
Particle p;
p.initialize();
std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz);
std::copy(dir, dir+3, p.coord[0].uvw);
p.coord[0].universe = model::root_universe;
p.r() = ll;
p.u() = u;
p.coord_[0].universe = model::root_universe;
// Open binary plot file for writing
std::ofstream of;
@ -891,9 +873,9 @@ void create_voxel(Plot pl)
voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
// move to center of voxels
ll[0] = ll[0] + vox[0] / 2.;
ll[1] = ll[1] + vox[1] / 2.;
ll[2] = ll[2] + vox[2] / 2.;
ll.x += vox[0] / 2.;
ll.y += vox[1] / 2.;
ll.z += vox[2] / 2.;
int data[pl.pixels_[1]][pl.pixels_[0]];
@ -910,16 +892,16 @@ void create_voxel(Plot pl)
// write to plot data
data[y][x] = id;
// advance particle in x direction
p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0];
p.r().x += vox[0];
}
// advance particle in y direction
p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1];
p.coord[0].xyz[0] = ll[0];
p.r().y += vox[1];
p.r().x = ll[0];
}
// advance particle in z direction
p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2];
p.coord[0].xyz[1] = ll[1];
p.coord[0].xyz[0] = ll[0];
p.r().z += vox[2];
p.r().y = ll[1];
p.r().x = ll[0];
// Write to HDF5 dataset
voxel_write_slice(z, dspace, dset, memspace, &(data[0]));
}

6
src/position.cpp
View file

@ -78,4 +78,10 @@ Position::operator/=(double v)
return *this;
}
Position
Position::operator-() const
{
return {-x, -y, -z};
}
} // namespace openmc

2
src/reaction.cpp
View file

@ -72,7 +72,7 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
// mark fission reactions so that we avoid the angle sampling.
if (is_fission(mt_)) {
for (auto& p : products_) {
if (p.particle_ == ParticleType::neutron) {
if (p.particle_ == Particle::Type::neutron) {
for (auto& d : p.distribution_) {
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) d_->fission() = true;

12
src/reaction_product.cpp
View file

@ -22,9 +22,9 @@ ReactionProduct::ReactionProduct(hid_t group)
std::string temp;
read_attribute(group, "particle", temp);
if (temp == "neutron") {
particle_ = ParticleType::neutron;
particle_ = Particle::Type::neutron;
} else if (temp == "photon") {
particle_ = ParticleType::photon;
particle_ = Particle::Type::photon;
}
// Read emission mode and decay rate
@ -69,13 +69,13 @@ ReactionProduct::ReactionProduct(hid_t group)
// Determine distribution type and read data
read_attribute(dgroup, "type", temp);
if (temp == "uncorrelated") {
distribution_.emplace_back(new UncorrelatedAngleEnergy{dgroup});
distribution_.push_back(std::make_unique<UncorrelatedAngleEnergy>(dgroup));
} else if (temp == "correlated") {
distribution_.emplace_back(new CorrelatedAngleEnergy{dgroup});
distribution_.push_back(std::make_unique<CorrelatedAngleEnergy>(dgroup));
} else if (temp == "nbody") {
distribution_.emplace_back(new NBodyPhaseSpace{dgroup});
distribution_.push_back(std::make_unique<NBodyPhaseSpace>(dgroup));
} else if (temp == "kalbach-mann") {
distribution_.emplace_back(new KalbachMann{dgroup});
distribution_.push_back(std::make_unique<KalbachMann>(dgroup));
}
close_group(dgroup);

4
src/settings.cpp
View file

@ -509,7 +509,7 @@ void read_settings_xml()
// Read entropy mesh from <entropy>
auto node_entropy = root.child("entropy");
model::meshes.emplace_back(new RegularMesh{node_entropy});
model::meshes.push_back(std::make_unique<RegularMesh>(node_entropy));
// Set entropy mesh index
index_entropy_mesh = model::meshes.size() - 1;
@ -555,7 +555,7 @@ void read_settings_xml()
// Read entropy mesh from <entropy>
auto node_ufs = root.child("uniform_fs");
model::meshes.emplace_back(new RegularMesh{node_ufs});
model::meshes.push_back(std::make_unique<RegularMesh>(node_ufs));
// Set entropy mesh index
index_ufs_mesh = model::meshes.size() - 1;

14
src/simulation.cpp
View file

@ -474,29 +474,29 @@ void initialize_history(Particle* p, int64_t index_source)
p->from_source(&simulation::source_bank[index_source - 1]);
// set identifier for particle
p->id = simulation::work_index[mpi::rank] + index_source;
p->id_ = simulation::work_index[mpi::rank] + index_source;
// set random number seed
int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
* settings::n_particles + p->id;
* settings::n_particles + p->id_;
set_particle_seed(particle_seed);
// set particle trace
simulation::trace = false;
if (simulation::current_batch == settings::trace_batch &&
simulation::current_gen == settings::trace_gen &&
p->id == settings::trace_particle) simulation::trace = true;
p->id_ == settings::trace_particle) simulation::trace = true;
// Set particle track.
p->write_track = false;
p->write_track_ = false;
if (settings::write_all_tracks) {
p->write_track = true;
p->write_track_ = true;
} else if (settings::track_identifiers.size() > 0) {
for (const auto& t : settings::track_identifiers) {
if (simulation::current_batch == t[0] &&
simulation::current_gen == t[1] &&
p->id == t[2]) {
p->write_track = true;
p->id_ == t[2]) {
p->write_track_ = true;
break;
}
}

29
src/source.cpp
View file

@ -47,9 +47,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
if (check_for_node(node, "particle")) {
auto temp_str = get_node_value(node, "particle", true, true);
if (temp_str == "neutron") {
particle_ = ParticleType::neutron;
particle_ = Particle::Type::neutron;
} else if (temp_str == "photon") {
particle_ = ParticleType::photon;
particle_ = Particle::Type::photon;
settings::photon_transport = true;
} else {
fatal_error(std::string("Unknown source particle type: ") + temp_str);
@ -140,9 +140,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
}
Bank SourceDistribution::sample() const
Particle::Bank SourceDistribution::sample() const
{
Bank site;
Particle::Bank site;
// Set weight to one by default
site.wgt = 1.0;
@ -153,17 +153,15 @@ Bank SourceDistribution::sample() const
static int n_accept = 0;
while (!found) {
// Set particle type
site.particle = static_cast<int>(particle_);
site.particle = particle_;
// Sample spatial distribution
Position r = space_->sample();
site.xyz[0] = r.x;
site.xyz[1] = r.y;
site.xyz[2] = r.z;
site.r = space_->sample();
double xyz[] {site.r.x, site.r.y, site.r.z};
// Now search to see if location exists in geometry
int32_t cell_index, instance;
int err = openmc_find_cell(site.xyz, &cell_index, &instance);
int err = openmc_find_cell(xyz, &cell_index, &instance);
found = (err != OPENMC_E_GEOMETRY);
// Check if spatial site is in fissionable material
@ -172,7 +170,7 @@ Bank SourceDistribution::sample() const
if (space_box) {
if (space_box->only_fissionable()) {
// Determine material
auto c = model::cells[cell_index];
const auto& c = model::cells[cell_index];
auto mat_index = c->material_.size() == 1
? c->material_[0] : c->material_[instance];
@ -200,10 +198,7 @@ Bank SourceDistribution::sample() const
++n_accept;
// Sample angle
Direction u = angle_->sample();
site.uvw[0] = u.x;
site.uvw[1] = u.y;
site.uvw[2] = u.z;
site.u = angle_->sample();
// Check for monoenergetic source above maximum particle energy
auto p = static_cast<int>(particle_);
@ -290,7 +285,7 @@ void initialize_source()
}
}
Bank sample_external_source()
Particle::Bank sample_external_source()
{
// Set the random number generator to the source stream.
prn_set_stream(STREAM_SOURCE);
@ -312,7 +307,7 @@ Bank sample_external_source()
}
// Sample source site from i-th source distribution
Bank site {model::external_sources[i].sample()};
Particle::Bank site {model::external_sources[i].sample()};
// If running in MG, convert site % E to group
if (!settings::run_CE) {

24
src/state_point.cpp
View file

@ -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<Bank> temp_source {simulation::source_bank.begin(),
std::vector<Particle::Bank> temp_source {simulation::source_bank.begin(),
simulation::source_bank.begin() + simulation::work};
#endif

8
src/summary.cpp
View file

@ -102,19 +102,19 @@ void write_geometry(hid_t file)
write_attribute(geom_group, "n_lattices", model::lattices.size());
auto cells_group = create_group(geom_group, "cells");
for (Cell* c : model::cells) c->to_hdf5(cells_group);
for (const auto& c : model::cells) c->to_hdf5(cells_group);
close_group(cells_group);
auto surfaces_group = create_group(geom_group, "surfaces");
for (Surface* surf : model::surfaces) surf->to_hdf5(surfaces_group);
for (const auto& surf : model::surfaces) surf->to_hdf5(surfaces_group);
close_group(surfaces_group);
auto universes_group = create_group(geom_group, "universes");
for (Universe* u : model::universes) u->to_hdf5(universes_group);
for (const auto& u : model::universes) u->to_hdf5(universes_group);
close_group(universes_group);
auto lattices_group = create_group(geom_group, "lattices");
for (Lattice* lat : model::lattices) lat->to_hdf5(lattices_group);
for (const auto& lat : model::lattices) lat->to_hdf5(lattices_group);
close_group(lattices_group);
close_group(geom_group);

39
src/surface.cpp
View file

@ -27,10 +27,8 @@ extern "C" const int BC_PERIODIC {3};
//==============================================================================
namespace model {
std::vector<Surface*> surfaces;
std::map<int, int> surface_map;
std::vector<std::unique_ptr<Surface>> surfaces;
std::unordered_map<int, int> surface_map;
} // namespace model
//==============================================================================
@ -1080,40 +1078,40 @@ void read_surfaces(pugi::xml_node node)
std::string surf_type = get_node_value(surf_node, "type", true, true);
if (surf_type == "x-plane") {
model::surfaces.push_back(new SurfaceXPlane(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceXPlane>(surf_node));
} else if (surf_type == "y-plane") {
model::surfaces.push_back(new SurfaceYPlane(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceYPlane>(surf_node));
} else if (surf_type == "z-plane") {
model::surfaces.push_back(new SurfaceZPlane(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceZPlane>(surf_node));
} else if (surf_type == "plane") {
model::surfaces.push_back(new SurfacePlane(surf_node));
model::surfaces.push_back(std::make_unique<SurfacePlane>(surf_node));
} else if (surf_type == "x-cylinder") {
model::surfaces.push_back(new SurfaceXCylinder(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceXCylinder>(surf_node));
} else if (surf_type == "y-cylinder") {
model::surfaces.push_back(new SurfaceYCylinder(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceYCylinder>(surf_node));
} else if (surf_type == "z-cylinder") {
model::surfaces.push_back(new SurfaceZCylinder(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceZCylinder>(surf_node));
} else if (surf_type == "sphere") {
model::surfaces.push_back(new SurfaceSphere(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceSphere>(surf_node));
} else if (surf_type == "x-cone") {
model::surfaces.push_back(new SurfaceXCone(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceXCone>(surf_node));
} else if (surf_type == "y-cone") {
model::surfaces.push_back(new SurfaceYCone(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceYCone>(surf_node));
} else if (surf_type == "z-cone") {
model::surfaces.push_back(new SurfaceZCone(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceZCone>(surf_node));
} else if (surf_type == "quadric") {
model::surfaces.push_back(new SurfaceQuadric(surf_node));
model::surfaces.push_back(std::make_unique<SurfaceQuadric>(surf_node));
} else {
std::stringstream err_msg;
@ -1143,8 +1141,8 @@ void read_surfaces(pugi::xml_node node)
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = model::surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
Surface* surf_base = model::surfaces[i_surf].get();
auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Make sure this surface inherits from PeriodicSurface.
if (!surf) {
@ -1188,8 +1186,8 @@ void read_surfaces(pugi::xml_node node)
for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
// Downcast to the PeriodicSurface type.
Surface* surf_base = model::surfaces[i_surf];
PeriodicSurface* surf = dynamic_cast<PeriodicSurface*>(surf_base);
Surface* surf_base = model::surfaces[i_surf].get();
auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
// Also try downcasting to the SurfacePlane type (which must be handled
// differently).
@ -1260,7 +1258,6 @@ void read_surfaces(pugi::xml_node node)
void free_memory_surfaces()
{
for (Surface* surf : model::surfaces) {delete surf;}
model::surfaces.clear();
model::surface_map.clear();
}

92
src/tallies/derivative.cpp
View file

@ -126,11 +126,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
if (score_bin == SCORE_FLUX) {
score *= flux_deriv;
return;
} else if (p->material == MATERIAL_VOID) {
} else if (p->material_ == MATERIAL_VOID) {
score *= flux_deriv;
return;
}
const Material& material {*model::materials[p->material]};
const Material& material {*model::materials[p->material_]};
if (material.id_ != deriv.diff_material) {
score *= flux_deriv;
return;
@ -190,7 +190,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
switch (tally.estimator_) {
case ESTIMATOR_ANALOG:
if (p->event_nuclide != deriv.diff_nuclide) {
if (p->event_nuclide_ != deriv.diff_nuclide) {
score *= flux_deriv;
return;
}
@ -320,10 +320,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
// Find the index of the event nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
if (material.nuclide_[i] == p->event_nuclide) break;
if (material.nuclide_[i] == p->event_nuclide_) break;
const auto& nuc {*data::nuclides[p->event_nuclide]};
if (!multipole_in_range(&nuc, p->last_E)) {
const auto& nuc {*data::nuclides[p->event_nuclide_]};
if (!multipole_in_range(&nuc, p->E_last_)) {
score *= flux_deriv;
break;
}
@ -331,10 +331,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
switch (score_bin) {
case SCORE_TOTAL:
if (simulation::micro_xs[p->event_nuclide].total) {
if (simulation::micro_xs[p->event_nuclide_].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
/ simulation::material_xs.total;
} else {
@ -343,11 +343,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_SCATTER:
if (simulation::micro_xs[p->event_nuclide].total
- simulation::micro_xs[p->event_nuclide].absorption) {
if (simulation::micro_xs[p->event_nuclide_].total
- simulation::micro_xs[p->event_nuclide_].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_s * material.atom_density_(i)
/ (simulation::material_xs.total
- simulation::material_xs.absorption);
@ -357,10 +357,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_ABSORPTION:
if (simulation::micro_xs[p->event_nuclide].absorption) {
if (simulation::micro_xs[p->event_nuclide_].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_a * material.atom_density_(i)
/ simulation::material_xs.absorption;
} else {
@ -369,10 +369,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_FISSION:
if (simulation::micro_xs[p->event_nuclide].fission) {
if (simulation::micro_xs[p->event_nuclide_].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_f * material.atom_density_(i)
/ simulation::material_xs.fission;
} else {
@ -381,12 +381,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
break;
case SCORE_NU_FISSION:
if (simulation::micro_xs[p->event_nuclide].fission) {
double nu = simulation::micro_xs[p->event_nuclide].nu_fission
/ simulation::micro_xs[p->event_nuclide].fission;
if (simulation::micro_xs[p->event_nuclide_].fission) {
double nu = simulation::micro_xs[p->event_nuclide_].nu_fission
/ simulation::micro_xs[p->event_nuclide_].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
/ simulation::material_xs.nu_fission;
} else {
@ -404,7 +404,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
case ESTIMATOR_COLLISION:
if (i_nuclide != -1) {
const auto& nuc {data::nuclides[i_nuclide]};
if (!multipole_in_range(nuc.get(), p->last_E)) {
if (!multipole_in_range(nuc.get(), p->E_last_)) {
score *= flux_deriv;
return;
}
@ -418,11 +418,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i);
}
}
@ -431,7 +431,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total;
} else {
@ -446,12 +446,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
if (multipole_in_range(&nuc, p->E_last_)
&& (simulation::micro_xs[i_nuc].total
- simulation::micro_xs[i_nuc].absorption)) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_s * material.atom_density_(i);
}
}
@ -462,7 +462,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
- simulation::micro_xs[i_nuclide].absorption);
} else {
@ -476,11 +476,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_a * material.atom_density_(i);
}
}
@ -489,7 +489,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_a / simulation::micro_xs[i_nuclide].absorption;
} else {
@ -503,11 +503,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += dsig_f * material.atom_density_(i);
}
}
@ -516,7 +516,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
} else {
@ -530,13 +530,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
if (multipole_in_range(&nuc, p->E_last_)
&& simulation::micro_xs[i_nuc].fission) {
double nu = simulation::micro_xs[i_nuc].nu_fission
/ simulation::micro_xs[i_nuc].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
cum_dsig += nu * dsig_f * material.atom_density_(i);
}
}
@ -545,7 +545,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
} else {
@ -570,8 +570,8 @@ void
score_track_derivative(const Particle* p, double distance)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material]};
if (p->material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;
@ -597,13 +597,13 @@ score_track_derivative(const Particle* p, double distance)
case DIFF_TEMPERATURE:
for (auto i = 0; i < material.nuclide_.size(); ++i) {
const auto& nuc {*data::nuclides[material.nuclide_[i]]};
if (multipole_in_range(&nuc, p->last_E)) {
if (multipole_in_range(&nuc, p->E_last_)) {
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
// (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_);
deriv.flux_deriv -= distance * (dsig_s + dsig_a)
* material.atom_density_(i);
}
@ -616,8 +616,8 @@ score_track_derivative(const Particle* p, double distance)
void score_collision_derivative(const Particle* p)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material]};
if (p->material_ == MATERIAL_VOID) return;
const Material& material {*model::materials[p->material_]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;
@ -632,7 +632,7 @@ void score_collision_derivative(const Particle* p)
break;
case DIFF_NUCLIDE_DENSITY:
if (p->event_nuclide != deriv.diff_nuclide) continue;
if (p->event_nuclide_ != deriv.diff_nuclide) continue;
// Find the index in this material for the diff_nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
@ -656,19 +656,19 @@ void score_collision_derivative(const Particle* p)
// Loop over the material's nuclides until we find the event nuclide.
for (auto i_nuc : material.nuclide_) {
const auto& nuc {*data::nuclides[i_nuc]};
if (i_nuc == p->event_nuclide && multipole_in_range(&nuc, p->last_E)) {
if (i_nuc == p->event_nuclide_ && multipole_in_range(&nuc, p->E_last_)) {
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
const auto& micro_xs {simulation::micro_xs[i_nuc]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
= nuc.multipole_->evaluate_deriv(p->E_last_, p->sqrtkT_);
deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
// Note that this is an approximation! The real scattering cross
// section is
// Sigma_s(E'->E, uvw'->uvw) = Sigma_s(E') * P(E'->E, uvw'->uvw).
// We are assuming that d_P(E'->E, uvw'->uvw) / d_T = 0 and only
// Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u).
// We are assuming that d_P(E'->E, u'->u) / d_T = 0 and only
// computing d_S(E') / d_T. Using this approximation in the vicinity
// of low-energy resonances causes errors (~2-5% for PWR pincell
// eigenvalue derivatives).

4
src/tallies/filter_azimuthal.cpp
View file

@ -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()) {

4
src/tallies/filter_cell.cpp
View file

@ -41,8 +41,8 @@ void
CellFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
for (int i = 0; i < p->n_coord; i++) {
auto search = map_.find(p->coord[i].cell);
for (int i = 0; i < p->n_coord_; i++) {
auto search = map_.find(p->coord_[i].cell);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

2
src/tallies/filter_cellborn.cpp
View file

@ -8,7 +8,7 @@ void
CellbornFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
auto search = map_.find(p->cell_born);
auto search = map_.find(p->cell_born_);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

4
src/tallies/filter_cellfrom.cpp
View file

@ -8,8 +8,8 @@ void
CellFromFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
for (int i = 0; i < p->last_n_coord; i++) {
auto search = map_.find(p->last_cell[i]);
for (int i = 0; i < p->n_coord_last_; i++) {
auto search = map_.find(p->cell_last_[i]);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

14
src/tallies/filter_distribcell.cpp
View file

@ -40,20 +40,20 @@ DistribcellFilter::get_all_bins(const Particle* p, int estimator,
{
int offset = 0;
auto distribcell_index = model::cells[cell_]->distribcell_index_;
for (int i = 0; i < p->n_coord; i++) {
auto& c {*model::cells[p->coord[i].cell]};
for (int i = 0; i < p->n_coord_; i++) {
auto& c {*model::cells[p->coord_[i].cell]};
if (c.type_ == FILL_UNIVERSE) {
offset += c.offset_[distribcell_index];
} else if (c.type_ == FILL_LATTICE) {
auto& lat {*model::lattices[p->coord[i+1].lattice-1]};
int i_xyz[3] {p->coord[i+1].lattice_x,
p->coord[i+1].lattice_y,
p->coord[i+1].lattice_z};
auto& lat {*model::lattices[p->coord_[i+1].lattice]};
int i_xyz[3] {p->coord_[i+1].lattice_x,
p->coord_[i+1].lattice_y,
p->coord_[i+1].lattice_z};
if (lat.are_valid_indices(i_xyz)) {
offset += lat.offset(distribcell_index, i_xyz);
}
}
if (cell_ == p->coord[i].cell) {
if (cell_ == p->coord_[i].cell) {
match.bins_.push_back(offset);
match.weights_.push_back(1.0);
return;

16
src/tallies/filter_energy.cpp
View file

@ -41,17 +41,17 @@ void
EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const
{
if (p->g != F90_NONE && matches_transport_groups_) {
if (p->g_ != F90_NONE && matches_transport_groups_) {
if (estimator == ESTIMATOR_TRACKLENGTH) {
match.bins_.push_back(data::num_energy_groups - p->g);
match.bins_.push_back(data::num_energy_groups - p->g_);
} else {
match.bins_.push_back(data::num_energy_groups - p->last_g);
match.bins_.push_back(data::num_energy_groups - p->g_last_);
}
match.weights_.push_back(1.0);
} else {
// Get the pre-collision energy of the particle.
auto E = p->last_E;
auto E = p->E_last_;
// Bin the energy.
if (E >= bins_.front() && E <= bins_.back()) {
@ -85,13 +85,13 @@ void
EnergyoutFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
if (p->g != F90_NONE && matches_transport_groups_) {
match.bins_.push_back(data::num_energy_groups - p->g);
if (p->g_ != F90_NONE && matches_transport_groups_) {
match.bins_.push_back(data::num_energy_groups - p->g_);
match.weights_.push_back(1.0);
} else {
if (p->E >= bins_.front() && p->E <= bins_.back()) {
auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E);
if (p->E_ >= bins_.front() && p->E_ <= bins_.back()) {
auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E_);
match.bins_.push_back(bin);
match.weights_.push_back(1.0);
}

6
src/tallies/filter_energyfunc.cpp
View file

@ -33,12 +33,12 @@ void
EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) {
if (p->E_last_ >= energy_.front() && p->E_last_ <= energy_.back()) {
// Search for the incoming energy bin.
auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E);
auto i = lower_bound_index(energy_.begin(), energy_.end(), p->E_last_);
// Compute the interpolation factor between the nearest bins.
double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]);
double f = (p->E_last_ - energy_[i]) / (energy_[i+1] - energy_[i]);
// Interpolate on the lin-lin grid.
match.bins_.push_back(0);

2
src/tallies/filter_legendre.cpp
View file

@ -19,7 +19,7 @@ LegendreFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
double wgt[n_bins_];
calc_pn_c(order_, p->mu, wgt);
calc_pn_c(order_, p->mu_, wgt);
for (int i = 0; i < n_bins_; i++) {
match.bins_.push_back(i);
match.weights_.push_back(wgt[i]);

2
src/tallies/filter_material.cpp
View file

@ -42,7 +42,7 @@ void
MaterialFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
auto search = map_.find(p->material);
auto search = map_.find(p->material_);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

2
src/tallies/filter_mesh.cpp
View file

@ -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);

4
src/tallies/filter_mu.cpp
View file

@ -38,8 +38,8 @@ void
MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const
{
if (p->mu >= bins_.front() && p->mu <= bins_.back()) {
auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu);
if (p->mu_ >= bins_.front() && p->mu_ <= bins_.back()) {
auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu_);
match.bins_.push_back(bin);
match.weights_.push_back(1.0);
}

25
src/tallies/filter_particle.cpp
View file

@ -7,8 +7,10 @@ namespace openmc {
void
ParticleFilter::from_xml(pugi::xml_node node)
{
particles_ = get_node_array<int>(node, "bins");
for (auto& p : particles_) --p;
auto particles = get_node_array<int>(node, "bins");
for (auto& p : particles) {
particles_.push_back(static_cast<Particle::Type>(p - 1));
}
n_bins_ = particles_.size();
}
@ -17,7 +19,7 @@ ParticleFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
for (auto i = 0; i < particles_.size(); i++) {
if (particles_[i] == p->type) {
if (particles_[i] == p->type_) {
match.bins_.push_back(i);
match.weights_.push_back(1.0);
}
@ -28,13 +30,26 @@ void
ParticleFilter::to_statepoint(hid_t filter_group) const
{
Filter::to_statepoint(filter_group);
write_dataset(filter_group, "bins", particles_);
std::vector<int> particles;
for (auto p : particles_) {
particles.push_back(static_cast<int>(p) + 1);
}
write_dataset(filter_group, "bins", particles);
}
std::string
ParticleFilter::text_label(int bin) const
{
return "Particle " + std::to_string(particles_[bin]);
switch (particles_[bin]) {
case Particle::Type::neutron:
return "Particle: neutron";
case Particle::Type::photon:
return "Particle: photon";
case Particle::Type::electron:
return "Particle: electron";
case Particle::Type::positron:
return "Particle: positron";
}
}
} // namespace openmc

4
src/tallies/filter_polar.cpp
View file

@ -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()) {

4
src/tallies/filter_sph_harm.cpp
View file

@ -37,14 +37,14 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator,
// Determine cosine term for scatter expansion if necessary
double wgt[order_ + 1];
if (cosine_ == SphericalHarmonicsCosine::scatter) {
calc_pn_c(order_, p->mu, wgt);
calc_pn_c(order_, p->mu_, wgt);
} else {
for (int i = 0; i < order_ + 1; i++) wgt[i] = 1;
}
// Find the Rn,m values
double rn[n_bins_];
calc_rn_c(order_, p->last_uvw, rn);
calc_rn(order_, p->u_last_, rn);
int j = 0;
for (int n = 0; n < order_ + 1; n++) {

6
src/tallies/filter_sptl_legendre.cpp
View file

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

4
src/tallies/filter_surface.cpp
View file

@ -41,10 +41,10 @@ void
SurfaceFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
auto search = map_.find(std::abs(p->surface)-1);
auto search = map_.find(std::abs(p->surface_)-1);
if (search != map_.end()) {
match.bins_.push_back(search->second);
if (p->surface < 0) {
if (p->surface_ < 0) {
match.weights_.push_back(-1.0);
} else {
match.weights_.push_back(1.0);

4
src/tallies/filter_universe.cpp
View file

@ -41,8 +41,8 @@ void
UniverseFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
for (int i = 0; i < p->n_coord; i++) {
auto search = map_.find(p->coord[i].universe);
for (int i = 0; i < p->n_coord_; i++) {
auto search = map_.find(p->coord_[i].universe);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

8
src/tallies/filter_zernike.cpp
View file

@ -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) {

12
src/tallies/tally.cpp
View file

@ -765,9 +765,9 @@ void read_tallies_xml()
} else {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (int p : pf->particles_) {
if (p == static_cast<int>(ParticleType::electron) ||
p == static_cast<int>(ParticleType::positron)) {
for (auto p : pf->particles_) {
if (p == Particle::Type::electron ||
p == Particle::Type::positron) {
t->estimator_ = ESTIMATOR_ANALOG;
}
}
@ -776,11 +776,11 @@ void read_tallies_xml()
if (particle_filter_index >= 0) {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (int p : pf->particles_) {
if (p != static_cast<int>(ParticleType::neutron)) {
for (auto p : pf->particles_) {
if (p != Particle::Type::neutron) {
warning("Particle filter other than NEUTRON used with photon "
"transport turned off. All tallies for particle type " +
std::to_string(p) + " will have no scores");
std::to_string(static_cast<int>(p)) + " will have no scores");
}
}
}

539
src/tallies/tally_scoring.cpp
View file

File diff suppressed because it is too large Load diff

4
src/track_output.cpp
View file

@ -40,14 +40,14 @@ void add_particle_track()
void write_particle_track(const Particle& p)
{
tracks.back().push_back({p.coord[0].xyz});
tracks.back().push_back(p.r());
}
void finalize_particle_track(const Particle& p)
{
std::stringstream filename;
filename << settings::path_output << "track_" << simulation::current_batch
<< '_' << simulation::current_gen << '_' << p.id << ".h5";
<< '_' << simulation::current_gen << '_' << p.id_ << ".h5";
// Determine number of coordinates for each particle
std::vector<int> n_coords;

29
src/volume_calc.cpp
View file

@ -96,7 +96,6 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
std::vector<std::vector<int>> indices(n);
std::vector<std::vector<int>> hits(n);
Particle p;
p.initialize();
prn_set_stream(STREAM_VOLUME);
@ -105,41 +104,37 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
for (int i = i_start; i < i_end; i++) {
set_particle_seed(i);
p.n_coord = 1;
p.n_coord_ = 1;
Position xi {prn(), prn(), prn()};
Position r {lower_left_ + xi*(upper_right_ - lower_left_)};
// TODO: assign directly when xyz is Position
std::copy(&r.x, &r.x + 3, p.coord[0].xyz);
p.coord[0].uvw[0] = 0.5;
p.coord[1].uvw[1] = 0.5;
p.coord[2].uvw[2] = 0.5;
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};
// If this location is not in the geometry at all, move on to next block
if (!find_cell(&p, false)) continue;
if (domain_type_ == FILTER_MATERIAL) {
if (p.material != MATERIAL_VOID) {
if (p.material_ != MATERIAL_VOID) {
for (int i_domain = 0; i_domain < n; i_domain++) {
if (model::materials[p.material]->id_ == domain_ids_[i_domain]) {
this->check_hit(p.material, indices[i_domain], hits[i_domain]);
if (model::materials[p.material_]->id_ == domain_ids_[i_domain]) {
this->check_hit(p.material_, indices[i_domain], hits[i_domain]);
break;
}
}
}
} else if (domain_type_ == FILTER_CELL) {
for (int level = 0; level < p.n_coord; ++level) {
for (int level = 0; level < p.n_coord_; ++level) {
for (int i_domain=0; i_domain < n; i_domain++) {
if (model::cells[p.coord[level].cell]->id_ == domain_ids_[i_domain]) {
this->check_hit(p.material, indices[i_domain], hits[i_domain]);
if (model::cells[p.coord_[level].cell]->id_ == domain_ids_[i_domain]) {
this->check_hit(p.material_, indices[i_domain], hits[i_domain]);
break;
}
}
}
} else if (domain_type_ == FILTER_UNIVERSE) {
for (int level = 0; level < p.n_coord; ++level) {
for (int level = 0; level < p.n_coord_; ++level) {
for (int i_domain = 0; i_domain < n; ++i_domain) {
if (model::universes[p.coord[level].universe]->id_ == domain_ids_[i_domain]) {
check_hit(p.material, indices[i_domain], hits[i_domain]);
if (model::universes[p.coord_[level].universe]->id_ == domain_ids_[i_domain]) {
check_hit(p.material_, indices[i_domain], hits[i_domain]);
break;
}
}

2
tests/regression_tests/sourcepoint_batch/test.py
View file

@ -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"

1
tests/unit_tests/test_capi.py
View file

@ -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):