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

Use Position and Direction in Particle class

Browse source
This commit is contained in:
Paul Romano 2019-03-01 14:46:21 -06:00
parent 477309c917
commit 368f89697d
45 changed files with 423 additions and 487 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);

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>

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

51
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,15 +33,19 @@ constexpr int MAX_LOST_PARTICLES {10};
// Maximum number of lost particles, relative to the total number of particles
constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
//==============================================================================
// Class declarations
//==============================================================================
struct LocalCoord {
Position r; //!< particle position
Direction u; //!< particle direction
int cell {-1};
int universe {-1};
int lattice {-1};
int lattice_x {-1};
int lattice_y {-1};
int lattice_z {-1};
double xyz[3]; //!< particle position
double uvw[3]; //!< particle direction
bool rotated {false}; //!< Is the level rotated?
//! clear data from a single coordinate level
@ -59,6 +63,16 @@ public:
neutron, photon, electron, positron
};
//! Saved ("banked") state of a particle
struct Bank {
Position r;
Direction u;
double E;
double wgt;
int delayed_group;
Type particle;
};
// Constructors
Particle();
@ -85,13 +99,13 @@ public:
bool alive_ {true}; //!< is particle alive?
// Other physical data
double last_xyz_current_[3]; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
double last_xyz_[3]; //!< previous coordinates
double last_uvw_[3]; //!< previous direction coordinates
double last_wgt_ {1.0}; //!< pre-collision particle weight
double absorb_wgt_ {0.0}; //!< weight absorbed for survival biasing
Position r_last_current_; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
Position r_last_; //!< previous coordinates
Direction u_last_; //!< previous direction coordinates
double last_wgt_ {1.0}; //!< pre-collision particle weight
double absorb_wgt_ {0.0}; //!< weight absorbed for survival biasing
// What event took place
bool fission_ {false}; //!< did particle cause implicit fission
@ -126,6 +140,18 @@ public:
int64_t n_secondary_ {};
Bank secondary_bank_[MAX_SECONDARY];
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
//! resets all coordinate levels for the particle
void clear();
@ -133,11 +159,10 @@ public:
//
//! stores the current phase space attributes of the particle in the
//! secondary bank and increments the number of sites in the secondary bank.
//! \param uvw Direction of the secondary particle
//! \param u Direction of the secondary particle
//! \param E Energy of the secondary particle in [eV]
//! \param type Particle type
//! \param run_CE Whether continuous-energy data is being used
void create_secondary(const double* uvw, double E, Type type, bool run_CE);
void create_secondary(Direction u, double E, Type type);
//! initialize from a source site
//

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

5
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,7 +39,7 @@ public:
//! Sample from the external source distribution
//! \return Sampled site
Bank sample() const;
Particle::Bank sample() const;
// Properties
double strength() const { return strength_; }
@ -62,7 +61,7 @@ extern "C" void initialize_source();
//! Sample a site from all external source distributions in proportion to their
//! source strength
//! \return Sampled source site
Bank sample_external_source();
Particle::Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();

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.");

2
src/bremsstrahlung.cpp
View file

@ -108,7 +108,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost)
if (w > settings::energy_cutoff[photon]) {
// Create secondary photon
p.create_secondary(p.coord_[0].uvw, w, Particle::Type::photon, true);
p.create_secondary(p.u(), w, Particle::Type::photon);
*E_lost += w;
}
}

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!");

89
src/geometry.cpp
View file

@ -44,7 +44,7 @@ check_cell_overlap(Particle* p)
// Loop through each cell on this level
for (auto index_cell : univ.cells_) {
Cell& c = *model::cells[index_cell];
if (c.contains(p->coord_[j].xyz, p->coord_[j].uvw, p->surface_)) {
if (c.contains(p->coord_[j].r, p->coord_[j].u, p->surface_)) {
if (index_cell != p->coord_[j].cell) {
std::stringstream err_msg;
err_msg << "Overlapping cells detected: " << c.id_ << ", "
@ -78,8 +78,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->coord_[p->n_coord_-1].xyz};
Direction u {p->coord_[p->n_coord_-1].uvw};
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord_[p->n_coord_-1].cell = i_cell;
@ -99,8 +99,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
if (model::cells[i_cell]->universe_ != i_universe) continue;
// Check if this cell contains the particle.
Position r {p->coord_[p->n_coord_-1].xyz};
Direction u {p->coord_[p->n_coord_-1].uvw};
Position r {p->r_local()};
Direction u {p->u_local()};
auto surf = p->surface_;
if (model::cells[i_cell]->contains(r, u, surf)) {
p->coord_[p->n_coord_-1].cell = i_cell;
@ -169,36 +169,28 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
p->coord_[p->n_coord_].universe = c.fill_;
// Set the position and direction.
for (int i = 0; i < 3; i++) {
p->coord_[p->n_coord_].xyz[i] = p->coord_[p->n_coord_-1].xyz[i];
p->coord_[p->n_coord_].uvw[i] = p->coord_[p->n_coord_-1].uvw[i];
}
p->coord_[p->n_coord_].r = p->coord_[p->n_coord_-1].r;
p->coord_[p->n_coord_].u = p->coord_[p->n_coord_-1].u;
// Apply translation.
p->coord_[p->n_coord_].xyz[0] -= c.translation_.x;
p->coord_[p->n_coord_].xyz[1] -= c.translation_.y;
p->coord_[p->n_coord_].xyz[2] -= c.translation_.z;
p->coord_[p->n_coord_].r -= c.translation_;
// Apply rotation.
if (!c.rotation_.empty()) {
auto x = p->coord_[p->n_coord_].xyz[0];
auto y = p->coord_[p->n_coord_].xyz[1];
auto z = p->coord_[p->n_coord_].xyz[2];
p->coord_[p->n_coord_].xyz[0] = x*c.rotation_[3] + y*c.rotation_[4]
+ z*c.rotation_[5];
p->coord_[p->n_coord_].xyz[1] = x*c.rotation_[6] + y*c.rotation_[7]
+ z*c.rotation_[8];
p->coord_[p->n_coord_].xyz[2] = x*c.rotation_[9] + y*c.rotation_[10]
+ z*c.rotation_[11];
auto u = p->coord_[p->n_coord_].uvw[0];
auto v = p->coord_[p->n_coord_].uvw[1];
auto w = p->coord_[p->n_coord_].uvw[2];
p->coord_[p->n_coord_].uvw[0] = u*c.rotation_[3] + v*c.rotation_[4]
+ w*c.rotation_[5];
p->coord_[p->n_coord_].uvw[1] = u*c.rotation_[6] + v*c.rotation_[7]
+ w*c.rotation_[8];
p->coord_[p->n_coord_].uvw[2] = u*c.rotation_[9] + v*c.rotation_[10]
+ w*c.rotation_[11];
Position r = p->coord_[p->n_coord_].r;
p->coord_[p->n_coord_].r.x = r.x*c.rotation_[3] + r.y*c.rotation_[4]
+ r.z*c.rotation_[5];
p->coord_[p->n_coord_].r.y = r.x*c.rotation_[6] + r.y*c.rotation_[7]
+ r.z*c.rotation_[8];
p->coord_[p->n_coord_].r.z = r.x*c.rotation_[9] + r.y*c.rotation_[10]
+ r.z*c.rotation_[11];
Direction u = p->coord_[p->n_coord_].u;
p->coord_[p->n_coord_].u.x = u.x*c.rotation_[3] + u.y*c.rotation_[4]
+ u.z*c.rotation_[5];
p->coord_[p->n_coord_].u.y = u.x*c.rotation_[6] + u.y*c.rotation_[7]
+ u.z*c.rotation_[8];
p->coord_[p->n_coord_].u.z = u.x*c.rotation_[9] + u.y*c.rotation_[10]
+ u.z*c.rotation_[11];
p->coord_[p->n_coord_].rotated = true;
}
@ -213,18 +205,12 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
Lattice& lat {*model::lattices[c.fill_]};
// Determine lattice indices.
Position r {p->coord_[p->n_coord_-1].xyz};
Direction u {p->coord_[p->n_coord_-1].uvw};
auto i_xyz = lat.get_indices(r, u);
auto i_xyz = lat.get_indices(p->r_local(), p->u_local());
// Store lower level coordinates.
r = lat.get_local_position(p->coord_[p->n_coord_-1].xyz, i_xyz);
p->coord_[p->n_coord_].xyz[0] = r.x;
p->coord_[p->n_coord_].xyz[1] = r.y;
p->coord_[p->n_coord_].xyz[2] = r.z;
p->coord_[p->n_coord_].uvw[0] = u.x;
p->coord_[p->n_coord_].uvw[1] = u.y;
p->coord_[p->n_coord_].uvw[2] = u.z;
Position r = lat.get_local_position(p->r_local(), i_xyz);
p->coord_[p->n_coord_].r = r;
p->coord_[p->n_coord_].u = p->u_local();
// Set lattice indices.
p->coord_[p->n_coord_].lattice = c.fill_;
@ -324,10 +310,7 @@ cross_lattice(Particle* p, int lattice_translation[3])
p->coord_[p->n_coord_-1].lattice_z};
// Set the new coordinate position.
auto r = lat.get_local_position(p->coord_[p->n_coord_-2].xyz, i_xyz);
p->coord_[p->n_coord_-1].xyz[0] = r.x;
p->coord_[p->n_coord_-1].xyz[1] = r.y;
p->coord_[p->n_coord_-1].xyz[2] = r.z;
p->r_local() = lat.get_local_position(p->coord_[p->n_coord_-2].r, i_xyz);
if (!lat.are_valid_indices(i_xyz)) {
// The particle is outside the lattice. Search for it from the base coords.
@ -377,8 +360,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
// Loop over each coordinate level.
for (int i = 0; i < p->n_coord_; i++) {
Position r {p->coord_[i].xyz};
Direction u {p->coord_[i].uvw};
Position r {p->coord_[i].r};
Direction u {p->coord_[i].u};
Cell& c {*model::cells[p->coord_[i].cell]};
// Find the oncoming surface in this cell and the distance to it.
@ -399,8 +382,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
lattice_distance = lat.distance(r, u, i_xyz);
break;
case LatticeType::hex:
Position r_hex {p->coord_[i-1].xyz[0], p->coord_[i-1].xyz[1],
p->coord_[i].xyz[2]};
Position r_hex {p->coord_[i-1].r.x, p->coord_[i-1].r.y,
p->coord_[i].r.z};
lattice_distance = lat.distance(r_hex, u, i_xyz);
break;
}
@ -466,15 +449,13 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
{
Particle p;
std::copy(xyz, xyz + 3, p.coord_[0].xyz);
p.coord_[0].uvw[0] = 0.0;
p.coord_[0].uvw[1] = 0.0;
p.coord_[0].uvw[2] = 1.0;
p.r() = Position{xyz};
p.u() = {0.0, 0.0, 1.0};
if (!find_cell(&p, false)) {
std::stringstream msg;
msg << "Could not find cell at position (" << xyz[0] << ", " << xyz[1]
<< ", " << xyz[2] << ").";
msg << "Could not find cell at position (" << p.r().x << ", " << p.r().y
<< ", " << p.r().z << ").";
set_errmsg(msg);
return OPENMC_E_GEOMETRY;
}

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

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

16
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;
@ -522,9 +522,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_;
@ -652,7 +652,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
@ -670,7 +670,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) {

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

4
src/mgxs_interface.cpp
View file

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

8
src/output.cpp
View file

@ -182,10 +182,10 @@ extern "C" void print_particle(Particle* p)
<< p->coord_[i].lattice_z << ")\n";
}
std::cout << " xyz = " << p->coord_[i].xyz[0] << " "
<< p->coord_[i].xyz[1] << " " << p->coord_[i].xyz[2] << "\n";
std::cout << " uvw = " << p->coord_[i].uvw[0] << " "
<< p->coord_[i].uvw[1] << " " << p->coord_[i].uvw[2] << "\n";
std::cout << " r = (" << p->coord_[i].r.x << ", "
<< p->coord_[i].r.y << ", " << p->coord_[i].r.z << ")\n";
std::cout << " u = (" << p->coord_[i].u.x << ", "
<< p->coord_[i].u.y << ", " << p->coord_[i].u.z << ")\n";
}
// Display miscellaneous info.

101
src/particle.cpp
View file

@ -66,19 +66,18 @@ Particle::clear()
}
void
Particle::create_secondary(const double* uvw, double E, Type type, bool run_CE)
Particle::create_secondary(Direction u, double E, Type type)
{
if (n_secondary_ == MAX_SECONDARY) {
fatal_error("Too many secondary particles created.");
}
int64_t n = n_secondary_;
secondary_bank_[n].particle = static_cast<int>(type);
secondary_bank_[n].particle = type;
secondary_bank_[n].wgt = wgt_;
std::copy(coord_[0].xyz, coord_[0].xyz + 3, secondary_bank_[n].xyz);
std::copy(uvw, uvw + 3, secondary_bank_[n].uvw);
secondary_bank_[n].E = E;
if (!run_CE) secondary_bank_[n].E = g_;
secondary_bank_[n].r = this->r();
secondary_bank_[n].u = u;
secondary_bank_[n].E = settings::run_CE ? E : g_;
++n_secondary_;
}
@ -95,14 +94,14 @@ Particle::from_source(const Bank* src)
fission_ = false;
// copy attributes from source bank site
type_ = static_cast<Particle::Type>(src->particle);
type_ = src->particle;
wgt_ = src->wgt;
last_wgt_ = src->wgt;
std::copy(src->xyz, src->xyz + 3, coord_[0].xyz);
std::copy(src->uvw, src->uvw + 3, coord_[0].uvw);
std::copy(src->xyz, src->xyz + 3, last_xyz_current_);
std::copy(src->xyz, src->xyz + 3, last_xyz_);
std::copy(src->uvw, src->uvw + 3, last_uvw_);
this->r() = src->r;
this->u() = src->u;
r_last_current_ = src->r;
r_last_ = src->r;
u_last_ = src->u;
if (settings::run_CE) {
E_ = src->E;
g_ = 0;
@ -153,8 +152,8 @@ Particle::transport()
// Store pre-collision particle properties
last_wgt_ = wgt_;
last_E_ = E_;
std::copy(coord_[0].uvw, coord_[0].uvw + 3, last_uvw_);
std::copy(coord_[0].xyz, coord_[0].xyz + 3, last_xyz_);
u_last_ = this->u();
r_last_ = this->r();
// If the cell hasn't been determined based on the particle's location,
// initiate a search for the current cell. This generally happens at the
@ -186,7 +185,7 @@ Particle::transport()
}
} else {
// Get the MG data
calculate_xs_c(material_, g_, sqrtkT_, coord_[n_coord_-1].uvw,
calculate_xs_c(material_, g_, sqrtkT_, this->u_local(),
simulation::material_xs.total, simulation::material_xs.absorption,
simulation::material_xs.nu_fission);
@ -225,10 +224,7 @@ Particle::transport()
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
// TODO: use Position
coord_[j].xyz[0] += distance * coord_[j].uvw[0];
coord_[j].xyz[1] += distance * coord_[j].uvw[1];
coord_[j].xyz[2] += distance * coord_[j].uvw[2];
coord_[j].r += distance * coord_[j].u;
}
// Score track-length tallies
@ -323,25 +319,25 @@ Particle::transport()
fission_ = false;
// Save coordinates for tallying purposes
std::copy(coord_[0].xyz, coord_[0].xyz + 3, last_xyz_current_);
r_last_current_ = this->r();
// Set last material to none since cross sections will need to be
// re-evaluated
last_material_ = C_NONE;
// Set all uvws to base level -- right now, after a collision, only the
// base level uvws are changed
// Set all directions to base level -- right now, after a collision, only
// the base level directions are changed
for (int j = 0; j < n_coord_ - 1; ++j) {
if (coord_[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord_[j].cell]->rotation_};
Direction u {coord_[j].uvw};
coord_[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord_[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord_[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z;
const auto& u {coord_[j].u};
coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
} else {
// Otherwise, copy this level's direction
std::copy(coord_[j].uvw, coord_[j].uvw + 3, coord_[j + 1].uvw);
coord_[j+1].u = coord_[j].u;
}
}
@ -403,10 +399,7 @@ Particle::cross_surface()
// TODO: Find a better solution to score surface currents than
// physically moving the particle forward slightly
// TODO: Use Position
coord_[0].xyz[0] += TINY_BIT * coord_[0].uvw[0];
coord_[0].xyz[1] += TINY_BIT * coord_[0].uvw[1];
coord_[0].xyz[2] += TINY_BIT * coord_[0].uvw[2];
this->r() += TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
}
@ -443,22 +436,17 @@ Particle::cross_surface()
if (!model::active_meshsurf_tallies.empty()) {
Position r {coord_[0].xyz};
coord_[0].xyz[0] -= TINY_BIT * coord_[0].uvw[0];
coord_[0].xyz[1] -= TINY_BIT * coord_[0].uvw[1];
coord_[0].xyz[2] -= TINY_BIT * coord_[0].uvw[2];
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
std::copy(&r.x, &r.x + 3, coord_[0].xyz);
this->r() = r;
}
// Reflect particle off surface
Direction u = surf->reflect(coord_[0].xyz, coord_[0].uvw);
Direction u = surf->reflect(this->r(), this->u());
// Make sure new particle direction is normalized
double norm = u.norm();
coord_[0].uvw[0] = u.x/norm;
coord_[0].uvw[1] = u.y/norm;
coord_[0].uvw[2] = u.z/norm;
this->u() = u / u.norm();
// Reassign particle's cell and surface
coord_[0].cell = last_cell_[last_n_coord_ - 1];
@ -476,9 +464,7 @@ Particle::cross_surface()
}
// Set previous coordinate going slightly past surface crossing
last_xyz_current_[0] = coord_[0].xyz[0] + TINY_BIT*coord_[0].uvw[0];
last_xyz_current_[1] = coord_[0].xyz[1] + TINY_BIT*coord_[0].uvw[1];
last_xyz_current_[2] = coord_[0].xyz[2] + TINY_BIT*coord_[0].uvw[2];
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
@ -502,12 +488,10 @@ Particle::cross_surface()
// particle to change -- artificially move the particle slightly back in
// case the surface crossing is coincident with a mesh boundary
if (!model::active_meshsurf_tallies.empty()) {
Position r {coord_[0].xyz};
coord_[0].xyz[0] -= TINY_BIT * coord_[0].uvw[0];
coord_[0].xyz[1] -= TINY_BIT * coord_[0].uvw[1];
coord_[0].xyz[2] -= TINY_BIT * coord_[0].uvw[2];
Position r {this->r()};
this->r() -= TINY_BIT * this->u();
score_surface_tally(this, model::active_meshsurf_tallies);
std::copy(&r.x, &r.x + 3, coord_[0].xyz);
this->r() = r;
}
// Get a pointer to the partner periodic surface
@ -516,11 +500,7 @@ Particle::cross_surface()
model::surfaces[surf_p->i_periodic_].get());
// Adjust the particle's location and direction.
Position r {coord_[0].xyz};
Direction u {coord_[0].uvw};
bool rotational = other->periodic_translate(surf_p, r, u);
std::copy(&r.x, &r.x + 3, coord_[0].xyz);
std::copy(&u.x, &u.x + 3, coord_[0].uvw);
bool rotational = other->periodic_translate(surf_p, this->r(), this->u());
// Reassign particle's surface
// TODO: off-by-one
@ -538,9 +518,7 @@ Particle::cross_surface()
}
// Set previous coordinate going slightly past surface crossing
last_xyz_current_[0] = coord_[0].xyz[0] + TINY_BIT * coord_[0].uvw[0];
last_xyz_current_[1] = coord_[0].xyz[1] + TINY_BIT * coord_[0].uvw[1];
last_xyz_current_[2] = coord_[0].xyz[2] + TINY_BIT * coord_[0].uvw[2];
r_last_current_ = this->r() + TINY_BIT*this->u();
// Diagnostic message
if (settings::verbosity >= 10 || simulation::trace) {
@ -588,9 +566,7 @@ Particle::cross_surface()
// forward a tiny bit it should fix the problem.
n_coord_ = 1;
coord_[0].xyz[0] += TINY_BIT * coord_[0].uvw[0];
coord_[0].xyz[1] += TINY_BIT * coord_[0].uvw[1];
coord_[0].xyz[2] += TINY_BIT * coord_[0].uvw[2];
this->r() += TINY_BIT * this->u();
// Couldn't find next cell anywhere! This probably means there is an actual
// undefined region in the geometry.
@ -673,9 +649,8 @@ Particle::write_restart() const
int64_t i = simulation::current_work;
write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
hsize_t dims[] {3};
write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false);
write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false);
write_dataset(file_id, "xyz", simulation::source_bank[i-1].r);
write_dataset(file_id, "uvw", simulation::source_bank[i-1].u);
// Close file
file_close(file_id);

13
src/particle_restart.cpp
View file

@ -45,11 +45,8 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
p.type_ = static_cast<Particle::Type>(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, "xyz", p.r());
read_dataset(file_id, "uvw", p.u());
// Set energy group and average energy in multi-group mode
if (!settings::run_CE) {
@ -59,9 +56,9 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
// Set particle last attributes
p.last_wgt_ = p.wgt_;
std::copy(p.coord_[0].xyz, p.coord_[0].xyz + 3, p.last_xyz_current_);
std::copy(p.coord_[0].xyz, p.coord_[0].xyz + 3, p.last_xyz_);
std::copy(p.coord_[0].uvw, p.coord_[0].uvw + 3, p.last_uvw_);
p.r_last_current_ = p.r();
p.r_last_ = p.r();
p.u_last_ = p.u();
p.last_E_ = p.E_;
p.last_g_ = p.g_;

22
src/photon.cpp
View file

@ -653,12 +653,12 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
if (shell.n_transitions == 0) {
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<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;
p.create_secondary(uvw.data(), E, Particle::Type::photon, true);
p.create_secondary(u, E, Particle::Type::photon);
return;
}
@ -678,10 +678,10 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Sample angle isotropically
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<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);
@ -690,7 +690,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Non-radiative transition -- Auger/Coster-Kronig effect
// Create auger electron
p.create_secondary(uvw.data(), E, Particle::Type::electron, true);
p.create_secondary(u, E, Particle::Type::electron);
// Fill hole left by emitted auger electron
int i_hole = shell_map_.at(secondary);
@ -700,7 +700,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Radiative transition -- get X-ray energy
// Create fluorescent photon
p.create_secondary(uvw.data(), E, Particle::Type::photon, true);
p.create_secondary(u, E, Particle::Type::photon);
}
// Fill hole created by electron transitioning to the photoelectron hole

132
src/physics.cpp
View file

@ -136,8 +136,8 @@ void sample_neutron_reaction(Particle* p)
}
void
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array,
int64_t* size_bank, int64_t bank_capacity)
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Particle::Bank* bank_array, int64_t* size_bank, int64_t bank_capacity)
{
// TODO: Heat generation from fission
@ -181,12 +181,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_
p->fission_ = true;
for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) {
// Bank source neutrons by copying the particle data
bank_array[i].xyz[0] = p->coord_[0].xyz[0];
bank_array[i].xyz[1] = p->coord_[0].xyz[1];
bank_array[i].xyz[2] = p->coord_[0].xyz[2];
bank_array[i].r = p->r();
// Set that the bank particle is a neutron
bank_array[i].particle = static_cast<int>(Particle::Type::neutron);
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
@ -244,7 +242,7 @@ void sample_photon_reaction(Particle* p)
prob += micro.coherent;
if (prob > cutoff) {
double mu = element.rayleigh_scatter(alpha);
rotate_angle_c(p->coord_[0].uvw, mu, nullptr);
p->u() = rotate_angle(p->u(), mu, nullptr);
p->event_mt_ = COHERENT;
return;
}
@ -270,10 +268,8 @@ void sample_photon_reaction(Particle* p)
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
double phi = 2.0*PI*prn();
double uvw[3];
std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_electron, &phi);
p->create_secondary(uvw, E_electron, Particle::Type::electron, true);
Direction u = rotate_angle(p->u(), mu_electron, &phi);
p->create_secondary(u, E_electron, Particle::Type::electron);
// TODO: Compton subshell data does not match atomic relaxation data
// Allow electrons to fill orbital and produce auger electrons
@ -285,7 +281,7 @@ void sample_photon_reaction(Particle* p)
phi += PI;
p->E_ = alpha_out*MASS_ELECTRON_EV;
rotate_angle_c(p->coord_[0].uvw, mu, &phi);
p->u() = rotate_angle(p->u(), mu, &phi);
p->event_mt_ = INCOHERENT;
return;
}
@ -326,13 +322,13 @@ void sample_photon_reaction(Particle* p)
}
double phi = 2.0*PI*prn();
std::array<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
p->create_secondary(uvw.data(), E_electron, Particle::Type::electron, true);
p->create_secondary(u, E_electron, Particle::Type::electron);
// Allow electrons to fill orbital and produce auger electrons
// and fluorescent photons
@ -355,15 +351,12 @@ void sample_photon_reaction(Particle* p)
&mu_electron, &mu_positron);
// Create secondary electron
double uvw[3];
std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_electron, nullptr);
p->create_secondary(uvw, E_electron, Particle::Type::electron, true);
Direction u = rotate_angle(p->u(), mu_electron, nullptr);
p->create_secondary(u, E_electron, Particle::Type::electron);
// Create secondary positron
std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu_positron, nullptr);
p->create_secondary(uvw, E_positron, Particle::Type::positron, true);
u = rotate_angle(p->u(), mu_positron, nullptr);
p->create_secondary(u, E_positron, Particle::Type::positron);
p->event_mt_ = PAIR_PROD;
p->alive_ = false;
@ -396,18 +389,14 @@ void sample_positron_reaction(Particle* p)
// Sample angle isotropically
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
std::array<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
p->create_secondary(uvw.data(), MASS_ELECTRON_EV, Particle::Type::photon, true);
uvw[0] = -uvw[0];
uvw[1] = -uvw[1];
uvw[2] = -uvw[2];
p->create_secondary(uvw.data(), MASS_ELECTRON_EV, Particle::Type::photon, true);
p->create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon);
p->create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon);
p->E_ = 0.0;
p->alive_ = false;
@ -585,7 +574,7 @@ void absorption(Particle* p, int i_nuclide)
void scatter(Particle* p, int i_nuclide)
{
// copy incoming direction
Direction u_old {p->coord_[0].uvw};
Direction u_old {p->u()};
// Get pointer to nuclide and grid index/interpolation factor
const auto& nuc {data::nuclides[i_nuclide]};
@ -613,8 +602,8 @@ void scatter(Particle* p, int i_nuclide)
double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT,
&p->E_, p->coord_[0].uvw, &p->mu_, &p->wgt_);
elastic_scatter(i_nuclide, *nuc->reactions_[0], kT,
p->E_, p->u(), p->mu_);
p->event_mt_ = ELASTIC;
sampled = true;
@ -625,7 +614,7 @@ void scatter(Particle* p, int i_nuclide)
// =======================================================================
// S(A,B) SCATTERING
sab_scatter(i_nuclide, micro.index_sab, &p->E_, p->coord_[0].uvw, &p->mu_);
sab_scatter(i_nuclide, micro.index_sab, p->E_, p->u(), p->mu_);
p->event_mt_ = ELASTIC;
sampled = true;
@ -673,39 +662,33 @@ void scatter(Particle* p, int i_nuclide)
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
Direction u_new;
u_new.x = mu;
u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p->mu_ = u_old.dot(u_new);
// Change direction of particle
p->coord_[0].uvw[0] = u_new.x;
p->coord_[0].uvw[1] = u_new.y;
p->coord_[0].uvw[2] = u_new.z;
p->u().x = mu;
p->u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
p->u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p->mu_ = u_old.dot(p->u());
}
}
}
void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
double* uvw, double* mu_lab, double* wgt)
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double& E,
Direction& u, double& mu_lab)
{
// get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide]};
double vel = std::sqrt(*E);
double vel = std::sqrt(E);
double awr = nuc->awr_;
// Neutron velocity in LAB
Direction u {uvw};
Direction v_n = vel*u;
// Sample velocity of target nucleus
Direction v_t {};
if (!simulation::micro_xs[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), *E, u, v_n,
simulation::micro_xs[i_nuclide].elastic, kT, wgt);
v_t = sample_target_velocity(nuc.get(), E, u, v_n,
simulation::micro_xs[i_nuclide].elastic, kT);
}
// Velocity of center-of-mass
@ -720,10 +703,10 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
// Sample scattering angle, checking if it is an ncorrelated angle-energy
// distribution
double mu_cm;
auto& d = rx->products_[0].distribution_[0];
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<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;
}
@ -739,26 +722,23 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
// Transform back to LAB frame
v_n += v_cm;
*E = v_n.dot(v_n);
vel = std::sqrt(*E);
E = v_n.dot(v_n);
vel = std::sqrt(E);
// compute cosine of scattering angle in LAB frame by taking dot product of
// neutron's pre- and post-collision angle
*mu_lab = u.dot(v_n) / vel;
mu_lab = u.dot(v_n) / vel;
// Set energy and direction of particle in LAB frame
u = v_n / vel;
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
// Because of floating-point roundoff, it may be possible for mu_lab to be
// outside of the range [-1,1). In these cases, we just set mu_lab to exactly
// -1 or 1
if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab);
if (std::abs(mu_lab) > 1.0) mu_lab = std::copysign(1.0, mu_lab);
}
void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu)
void sab_scatter(int i_nuclide, int i_sab, double& E, Direction& u, double& mu)
{
// Determine temperature index
const auto& micro {simulation::micro_xs[i_nuclide]};
@ -766,15 +746,15 @@ void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu)
// Sample energy and angle
double E_out;
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, *E, &E_out, mu);
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, E, &E_out, &mu);
// Set energy to outgoing, change direction of particle
*E = E_out;
rotate_angle_c(uvw, *mu, nullptr);
E = E_out;
u = rotate_angle(u, mu, nullptr);
}
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_neut, double xs_eff, double kT, double* wgt)
Direction v_neut, double xs_eff, double kT)
{
// check if nuclide is a resonant scatterer
ResScatMethod sampling_method;
@ -973,7 +953,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT)
return vt * rotate_angle(u, mu, nullptr);
}
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site)
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site)
{
// Sample cosine of angle -- fission neutrons are always emitted
// isotropically. Sometimes in ACE data, fission reactions actually have
@ -983,9 +963,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
// Sample azimuthal angle uniformly in [0,2*pi)
double phi = 2.0*PI*prn();
site->uvw[0] = mu;
site->uvw[1] = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->uvw[2] = std::sqrt(1.0 - mu*mu) * std::sin(phi);
site->u.x = mu;
site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
// Determine total nu, delayed nu, and delayed neutron fraction
const auto& nuc {data::nuclides[i_nuclide]};
@ -1098,14 +1078,14 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
p->mu_ = mu;
// change direction of particle
rotate_angle_c(p->coord_[0].uvw, mu, nullptr);
p->u() = rotate_angle(p->u(), mu, nullptr);
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
if (std::floor(yield) == yield) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
p->create_secondary(p->coord_[0].uvw, p->E_, Particle::Type::neutron, true);
p->create_secondary(p->u(), p->E_, Particle::Type::neutron);
}
} else {
// Otherwise, change weight of particle based on yield
@ -1135,12 +1115,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
rx->products_[i_product].sample(p->E_, E, mu);
// Sample the new direction
double uvw[3];
std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw);
rotate_angle_c(uvw, mu, nullptr);
Direction u = rotate_angle(p->u(), mu, nullptr);
// Create the secondary photon
p->create_secondary(uvw, E, Particle::Type::photon, true);
p->create_secondary(u, E, Particle::Type::photon);
}
}

18
src/physics_mg.cpp
View file

@ -92,7 +92,7 @@ scatter(Particle* p)
p->g_ = gout + 1;
// Rotate the angle
rotate_angle_c(p->coord_[0].uvw, p->mu_, nullptr);
p->u() = rotate_angle(p->u(), p->mu_, nullptr);
// Update energy value for downstream compatability (in tallying)
p->E_ = data::energy_bin_avg[gout];
@ -102,7 +102,7 @@ scatter(Particle* p)
}
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size)
{
// TODO: Heat generation from fission
@ -152,12 +152,10 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
for (size_t i = static_cast<size_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>(Particle::Type::neutron);
bank_array[i].particle = Particle::Type::neutron;
// Set the weight of the fission bank site
bank_array[i].wgt = 1. / weight;
@ -168,16 +166,16 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn();
bank_array[i].uvw[0] = mu;
bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi);
bank_array[i].u.x = mu;
bank_array[i].u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
bank_array[i].u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
// Sample secondary energy distribution for the fission reaction and set
// the energy in the fission bank
int dg;
int gout;
data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout);
bank_array[i].E = static_cast<double>(gout + 1);
bank_array[i].E = gout + 1;
bank_array[i].delayed_group = dg + 1;
// Set the delayed group on the particle as well

118
src/plot.cpp
View file

@ -106,48 +106,48 @@ void create_ppm(Plot pl)
data.resize({width, height});
int in_i, out_i;
double xyz[3];
Position r;
switch(pl.basis_) {
case PlotBasis::xy :
in_i = 0;
out_i = 1;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1] + pl.width_[1] / 2.;
xyz[2] = pl.origin_[2];
r.x = pl.origin_[0] - pl.width_[0] / 2.;
r.y = pl.origin_[1] + pl.width_[1] / 2.;
r.z = pl.origin_[2];
break;
case PlotBasis::xz :
in_i = 0;
out_i = 2;
xyz[0] = pl.origin_[0] - pl.width_[0] / 2.;
xyz[1] = pl.origin_[1];
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
r.x = pl.origin_[0] - pl.width_[0] / 2.;
r.y = pl.origin_[1];
r.z = pl.origin_[2] + pl.width_[1] / 2.;
break;
case PlotBasis::yz :
in_i = 1;
out_i = 2;
xyz[0] = pl.origin_[0];
xyz[1] = pl.origin_[1] - pl.width_[0] / 2.;
xyz[2] = pl.origin_[2] + pl.width_[1] / 2.;
r.x = pl.origin_[0];
r.y = pl.origin_[1] - pl.width_[0] / 2.;
r.z = pl.origin_[2] + pl.width_[1] / 2.;
break;
}
double dir[3] = {0.5, 0.5, 0.5};
Direction u {0.5, 0.5, 0.5};
#pragma omp parallel
{
Particle p;
std::copy(xyz, xyz+3, p.coord_[0].xyz);
std::copy(dir, dir+3, p.coord_[0].uvw);
p.r() = r;
p.u() = u;
p.coord_[0].universe = model::root_universe;
#pragma omp for
for (int y = 0; y < height; y++) {
p.coord_[0].xyz[out_i] = xyz[out_i] - out_pixel * y;
p.r()[out_i] = r[out_i] - out_pixel * y;
for (int x = 0; x < width; x++) {
// local variables
RGBColor rgb;
int id;
p.coord_[0].xyz[in_i] = xyz[in_i] + in_pixel * x;
p.r()[in_i] = r[in_i] + in_pixel * x;
position_rgb(p, pl, rgb, id);
data(x,y) = rgb;
}
@ -308,11 +308,9 @@ void
Plot::set_origin(pugi::xml_node plot_node)
{
// Copy plotting origin
std::vector<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 "
@ -328,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 "
@ -339,9 +337,7 @@ Plot::set_width(pugi::xml_node plot_node)
} else if (PlotType::voxel == type_) {
if (pl_width.size() == 3) {
pl_width = get_node_array<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 "
@ -745,36 +741,27 @@ void draw_mesh_lines(Plot pl, ImageData& data)
break;
}
double xyz_ll_plot[3], xyz_ur_plot[3];
xyz_ll_plot[0] = pl.origin_[0];
xyz_ll_plot[1] = pl.origin_[1];
xyz_ll_plot[2] = pl.origin_[2];
Position ll_plot {pl.origin_};
Position ur_plot {pl.origin_};
xyz_ur_plot[0] = pl.origin_[0];
xyz_ur_plot[1] = pl.origin_[1];
xyz_ur_plot[2] = pl.origin_[2];
ll_plot[outer] -= pl.width_[0] / 2.;
ll_plot[inner] -= pl.width_[1] / 2.;
ur_plot[outer] += pl.width_[0] / 2.;
ur_plot[inner] += pl.width_[1] / 2.;
xyz_ll_plot[outer] = pl.origin_[outer] - pl.width_[0] / 2.;
xyz_ll_plot[inner] = pl.origin_[inner] - pl.width_[1] / 2.;
xyz_ur_plot[outer] = pl.origin_[outer] + pl.width_[0] / 2.;
xyz_ur_plot[inner] = pl.origin_[inner] + pl.width_[1] / 2.;
int width[3];
width[0] = xyz_ur_plot[0] - xyz_ll_plot[0];
width[1] = xyz_ur_plot[1] - xyz_ll_plot[1];
width[2] = xyz_ur_plot[2] - xyz_ll_plot[2];
Position width = ur_plot - ll_plot;
auto& m = model::meshes[pl.index_meshlines_mesh_];
int ijk_ll[3], ijk_ur[3];
bool in_mesh;
m->get_indices(Position(xyz_ll_plot), &(ijk_ll[0]), &in_mesh);
m->get_indices(Position(xyz_ur_plot), &(ijk_ur[0]), &in_mesh);
m->get_indices(ll_plot, &(ijk_ll[0]), &in_mesh);
m->get_indices(ur_plot, &(ijk_ur[0]), &in_mesh);
// Fortran/C++ index correction
ijk_ur[0]++; ijk_ur[1]++; ijk_ur[2]++;
double xyz_ll[3], xyz_ur[3];
Position r_ll, r_ur;
// sweep through all meshbins on this plane and draw borders
for (int i = ijk_ll[outer]; i <= ijk_ur[outer]; i++) {
for (int j = ijk_ll[inner]; j <= ijk_ur[inner]; j++) {
@ -782,20 +769,20 @@ void draw_mesh_lines(Plot pl, ImageData& data)
if (i > 0 && i <= m->shape_[outer] && j >0 && j <= m->shape_[inner] ) {
int outrange[3], inrange[3];
// get xyz's of lower left and upper right of this mesh cell
xyz_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
xyz_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
xyz_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
xyz_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
r_ll[outer] = m->lower_left_[outer] + m->width_[outer] * (i - 1);
r_ll[inner] = m->lower_left_[inner] + m->width_[inner] * (j - 1);
r_ur[outer] = m->lower_left_[outer] + m->width_[outer] * i;
r_ur[inner] = m->lower_left_[inner] + m->width_[inner] * j;
// map the xyz ranges to pixel ranges
double frac = (xyz_ll[outer] - xyz_ll_plot[outer]) / width[outer];
double frac = (r_ll[outer] - ll_plot[outer]) / width[outer];
outrange[0] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[outer] - xyz_ll_plot[outer]) / width[outer];
frac = (r_ur[outer] - ll_plot[outer]) / width[outer];
outrange[1] = int(frac * double(pl.pixels_[0]));
frac = (xyz_ur[inner] - xyz_ll_plot[inner]) / width[inner];
frac = (r_ur[inner] - ll_plot[inner]) / width[inner];
inrange[0] = int((1. - frac) * (double)pl.pixels_[1]);
frac = (xyz_ll[inner] - xyz_ll_plot[inner]) / width[inner];
frac = (r_ll[inner] - ll_plot[inner]) / width[inner];
inrange[1] = int((1. - frac) * (double)pl.pixels_[1]);
// draw lines
@ -845,16 +832,13 @@ void create_voxel(Plot pl)
vox[2] = pl.width_[2]/(double)pl.pixels_[2];
// initial particle position
std::array<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;
std::copy(ll.begin(), ll.begin()+ll.size(), p.coord_[0].xyz);
std::copy(dir, dir+3, p.coord_[0].uvw);
p.r() = ll;
p.u() = u;
p.coord_[0].universe = model::root_universe;
// Open binary plot file for writing
@ -889,9 +873,9 @@ void create_voxel(Plot pl)
voxel_init(file_id, &(dims[0]), &dspace, &dset, &memspace);
// move to center of voxels
ll[0] = ll[0] + vox[0] / 2.;
ll[1] = ll[1] + vox[1] / 2.;
ll[2] = ll[2] + vox[2] / 2.;
ll.x += vox[0] / 2.;
ll.y += vox[1] / 2.;
ll.z += vox[2] / 2.;
int data[pl.pixels_[1]][pl.pixels_[0]];
@ -908,16 +892,16 @@ void create_voxel(Plot pl)
// write to plot data
data[y][x] = id;
// advance particle in x direction
p.coord_[0].xyz[0] = p.coord_[0].xyz[0] + vox[0];
p.r().x += vox[0];
}
// advance particle in y direction
p.coord_[0].xyz[1] = p.coord_[0].xyz[1] + vox[1];
p.coord_[0].xyz[0] = ll[0];
p.r().y += vox[1];
p.r().x = ll[0];
}
// advance particle in z direction
p.coord_[0].xyz[2] = p.coord_[0].xyz[2] + vox[2];
p.coord_[0].xyz[1] = ll[1];
p.coord_[0].xyz[0] = ll[0];
p.r().z += vox[2];
p.r().y = ll[1];
p.r().x = ll[0];
// Write to HDF5 dataset
voxel_write_slice(z, dspace, dset, memspace, &(data[0]));
}

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

22
src/source.cpp
View file

@ -140,9 +140,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
}
Bank SourceDistribution::sample() const
Particle::Bank SourceDistribution::sample() const
{
Bank site;
Particle::Bank site;
// Set weight to one by default
site.wgt = 1.0;
@ -153,17 +153,14 @@ Bank SourceDistribution::sample() const
static int n_accept = 0;
while (!found) {
// Set particle type
site.particle = static_cast<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();
// Now search to see if location exists in geometry
int32_t cell_index, instance;
int err = openmc_find_cell(site.xyz, &cell_index, &instance);
int err = openmc_find_cell(&site.r.x, &cell_index, &instance);
found = (err != OPENMC_E_GEOMETRY);
// Check if spatial site is in fissionable material
@ -200,10 +197,7 @@ Bank SourceDistribution::sample() const
++n_accept;
// Sample angle
Direction u = angle_->sample();
site.uvw[0] = u.x;
site.uvw[1] = u.y;
site.uvw[2] = u.z;
site.u = angle_->sample();
// Check for monoenergetic source above maximum particle energy
auto p = static_cast<int>(particle_);
@ -290,7 +284,7 @@ void initialize_source()
}
}
Bank sample_external_source()
Particle::Bank sample_external_source()
{
// Set the random number generator to the source stream.
prn_set_stream(STREAM_SOURCE);
@ -312,7 +306,7 @@ Bank sample_external_source()
}
// Sample source site from i-th source distribution
Bank site {model::external_sources[i].sample()};
Particle::Bank site {model::external_sources[i].sample()};
// If running in MG, convert site % E to group
if (!settings::run_CE) {

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

4
src/tallies/derivative.cpp
View file

@ -667,8 +667,8 @@ void score_collision_derivative(const Particle* p)
deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
// Note that this is an approximation! The real scattering cross
// section is
// Sigma_s(E'->E, uvw'->uvw) = Sigma_s(E') * P(E'->E, uvw'->uvw).
// We are assuming that d_P(E'->E, uvw'->uvw) / d_T = 0 and only
// Sigma_s(E'->E, u'->u) = Sigma_s(E') * P(E'->E, u'->u).
// We are assuming that d_P(E'->E, u'->u) / d_T = 0 and only
// computing d_S(E') / d_T. Using this approximation in the vicinity
// of low-energy resonances causes errors (~2-5% for PWR pincell
// eigenvalue derivatives).

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

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

2
src/tallies/filter_sph_harm.cpp
View file

@ -44,7 +44,7 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator,
// Find the Rn,m values
double rn[n_bins_];
calc_rn_c(order_, p->last_uvw_, rn);
calc_rn(order_, p->u_last_, rn);
int j = 0;
for (int n = 0; n < order_ + 1; n++) {

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

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

17
src/tallies/tally_scoring.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/tallies/tally_scoring.h"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
@ -1223,7 +1224,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
auto& tally {*model::tallies[i_tally]};
// Set the direction and group to use with get_xs
const double* p_uvw;
Direction p_u;
int p_g;
if (tally.estimator_ == ESTIMATOR_ANALOG
|| tally.estimator_ == ESTIMATOR_COLLISION) {
@ -1233,40 +1234,40 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Then we either are alive and had a scatter (and so g changed),
// or are dead and g did not change
if (p->alive_) {
p_uvw = p->last_uvw_;
p_u = p->u_last_;
p_g = p->last_g_;
} else {
p_uvw = p->coord_[p->n_coord_-1].uvw;
p_u = p->u_local();
p_g = p->g_;
}
} else if (p->event_ == EVENT_SCATTER) {
// Then the energy group has been changed by the scattering routine
// meaning gin is now in p % last_g
p_uvw = p->last_uvw_;
p_u = p->u_last_;
p_g = p->last_g_;
} else {
// No scatter, no change in g.
p_uvw = p->coord_[p->n_coord_-1].uvw;
p_u = p->u_local();
p_g = p->g_;
}
} else {
// No actual collision so g has not changed.
p_uvw = p->coord_[p->n_coord_-1].uvw;
p_u = p->u_local();
p_g = p->g_;
}
// To significantly reduce de-referencing, point matxs to the macroscopic
// Mgxs for the material of interest
data::macro_xs[p->material_].set_angle_index(p_uvw);
data::macro_xs[p->material_].set_angle_index(p_u);
// Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide >= 0) {
// And since we haven't calculated this temperature index yet, do so now
data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT_);
data::nuclides_MG[i_nuclide].set_angle_index(p_uvw);
data::nuclides_MG[i_nuclide].set_angle_index(p_u);
}
for (auto i = 0; i < tally.scores_.size(); ++i) {

2
src/track_output.cpp
View file

@ -40,7 +40,7 @@ void add_particle_track()
void write_particle_track(const Particle& p)
{
tracks.back().push_back({p.coord_[0].xyz});
tracks.back().push_back(p.r());
}
void finalize_particle_track(const Particle& p)

8
src/volume_calc.cpp
View file

@ -106,12 +106,8 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
p.n_coord_ = 1;
Position xi {prn(), prn(), prn()};
Position r {lower_left_ + xi*(upper_right_ - lower_left_)};
// TODO: assign directly when xyz is Position
std::copy(&r.x, &r.x + 3, p.coord_[0].xyz);
p.coord_[0].uvw[0] = 0.5;
p.coord_[1].uvw[1] = 0.5;
p.coord_[2].uvw[2] = 0.5;
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};
// If this location is not in the geometry at all, move on to next block
if (!find_cell(&p, false)) continue;

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