#include "openmc/particle.h"

#include <algorithm> // copy, min

#include <cmath> // log, abs

#include <fmt/core.h>

#include "openmc/bank.h"

#include "openmc/capi.h"

#include "openmc/cell.h"

#include "openmc/collision_track.h"

#include "openmc/constants.h"

#include "openmc/dagmc.h"

#include "openmc/error.h"

#include "openmc/geometry.h"

#include "openmc/hdf5_interface.h"

#include "openmc/lattice.h"

#include "openmc/material.h"

#include "openmc/message_passing.h"

#include "openmc/mgxs_interface.h"

#include "openmc/nuclide.h"

#include "openmc/particle_data.h"

#include "openmc/photon.h"

#include "openmc/physics.h"

#include "openmc/physics_mg.h"

#include "openmc/random_lcg.h"

#include "openmc/settings.h"

#include "openmc/simulation.h"

#include "openmc/source.h"

#include "openmc/surface.h"

#include "openmc/tallies/derivative.h"

#include "openmc/tallies/tally.h"

#include "openmc/tallies/tally_scoring.h"

#include "openmc/track_output.h"

#include "openmc/weight_windows.h"

#ifdef OPENMC_DAGMC_ENABLED

#include "DagMC.hpp"

#endif

namespace openmc {

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

// Particle implementation

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

double Particle::speed(double E) const

{

// Determine mass in eV/c^2

double mass = this->mass();

// Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:

return C_LIGHT * std::sqrt(E * (E + 2 * mass)) / (E + mass);

}

double Particle::speed() const

{

if (settings::run_CE) {

return speed(this->E());

} else {

auto mat = this->material();

if (mat == MATERIAL_VOID)

return 1.0 / data::mg.default_inverse_velocity_[this->g()];

auto& macro_xs = data::mg.macro_xs_[mat];

int macro_t = this->mg_xs_cache().t;

int macro_a = macro_xs.get_angle_index(this->u());

return 1.0 / macro_xs.get_xs(

MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a);

}

}

double Particle::mass() const

{

switch (type().pdg_number()) {

case PDG_NEUTRON:

return MASS_NEUTRON_EV;

case PDG_ELECTRON:

case PDG_POSITRON:

return MASS_ELECTRON_EV;

case PDG_PHOTON:

return 0.0;

default:

return this->type().mass() * AMU_EV;

}

}

bool Particle::create_secondary(

double wgt, Direction u, double E, ParticleType type)

{

// If energy is below cutoff for this particle, don't create secondary

// particle

int idx = type.transport_index();

if (idx == C_NONE) {

return false;

}

if (E < settings::energy_cutoff[idx]) {

return false;

}

// Increment number of secondaries created (for ParticleProductionFilter)

n_secondaries()++;

SourceSite bank;

bank.particle = type;

bank.wgt = wgt;

bank.r = r();

bank.u = u;

bank.E = settings::run_CE ? E : g();

bank.time = time();

bank_second_E() += bank.E;

bank.parent_id = current_work();

if (settings::use_shared_secondary_bank) {

bank.progeny_id = n_progeny()++;

}

bank.wgt_born = wgt_born();

bank.wgt_ww_born = wgt_ww_born();

bank.n_split = n_split();

local_secondary_bank().emplace_back(bank);

return true;

}

void Particle::split(double wgt)

{

SourceSite bank;

bank.particle = type();

bank.wgt = wgt;

bank.r = r();

bank.u = u();

bank.E = settings::run_CE ? E() : g();

bank.time = time();

// Convert signed index to a signed surface ID

if (surface() == SURFACE_NONE) {

bank.surf_id = SURFACE_NONE;

} else {

int surf_id = model::surfaces[surface_index()]->id_;

bank.surf_id = (surface() > 0) ? surf_id : -surf_id;

}

bank.wgt_born = wgt_born();

bank.wgt_ww_born = wgt_ww_born();

bank.n_split = n_split();

bank.n_collision = n_collision();

bank.parent_id = current_work();

if (settings::use_shared_secondary_bank) {

bank.progeny_id = n_progeny()++;

}

local_secondary_bank().emplace_back(bank);

}

void Particle::from_source(const SourceSite* src)

{

// Reset some attributes

clear();

surface() = SURFACE_NONE;

cell_born() = C_NONE;

material() = C_NONE;

n_collision() = src->n_collision;

fission() = false;

zero_flux_derivs();

lifetime() = 0.0;

#ifdef OPENMC_DAGMC_ENABLED

history().reset();

#endif

// Copy attributes from source bank site

type() = src->particle;

wgt() = src->wgt;

wgt_last() = src->wgt;

r() = src->r;

u() = src->u;

r_born() = src->r;

r_last_current() = src->r;

r_last() = src->r;

u_last() = src->u;

if (settings::run_CE) {

E() = src->E;

g() = 0;

} else {

g() = static_cast<int>(src->E);

g_last() = static_cast<int>(src->E);

E() = data::mg.energy_bin_avg_[g()];

}

E_last() = E();

time() = src->time;

time_last() = src->time;

parent_nuclide() = src->parent_nuclide;

delayed_group() = src->delayed_group;

// Convert signed surface ID to signed index

if (src->surf_id != SURFACE_NONE) {

auto it = model::surface_map.find(std::abs(src->surf_id));

if (it != model::surface_map.end()) {

int index_plus_one = it->second + 1;

surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;

}

}

wgt_born() = src->wgt_born;

wgt_ww_born() = src->wgt_ww_born;

n_split() = src->n_split;

}

void Particle::event_calculate_xs()

{

// Set the random number stream

stream() = STREAM_TRACKING;

// Store pre-collision particle properties

wgt_last() = wgt();

E_last() = E();

u_last() = u();

r_last() = r();

time_last() = time();

// Reset event variables

event() = TallyEvent::KILL;

event_nuclide() = NUCLIDE_NONE;

event_mt() = REACTION_NONE;

// 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 (lowest_coord().cell() == C_NONE) {

if (!exhaustive_find_cell(*this)) {

mark_as_lost(

"Could not find the cell containing particle " + std::to_string(id()));

return;

}

// Set birth cell attribute

if (cell_born() == C_NONE)

cell_born() = lowest_coord().cell();

// Initialize last cells from current cell

for (int j = 0; j < n_coord(); ++j) {

cell_last(j) = coord(j).cell();

}

n_coord_last() = n_coord();

}

// Write particle track.

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 (settings::run_CE) {

if (material() != material_last() || sqrtkT() != sqrtkT_last() ||

density_mult() != density_mult_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);

}

} else {

// Get the MG data; unlike the CE case above, we have to re-calculate

// cross sections for every collision since the cross sections may

// be angle-dependent

data::mg.macro_xs_[material()].calculate_xs(*this);

// Update the particle's group while we know we are multi-group

g_last() = g();

}

} else {

macro_xs().total = 0.0;

macro_xs().absorption = 0.0;

macro_xs().fission = 0.0;

macro_xs().nu_fission = 0.0;

}

}

void Particle::event_advance()

{

// Find the distance to the nearest boundary

boundary() = distance_to_boundary(*this);

// Sample a distance to collision

if (type() == ParticleType::electron() ||

type() == ParticleType::positron()) {

collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;

} else if (macro_xs().total == 0.0) {

collision_distance() = INFINITY;

} else {

collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;

}

double speed = this->speed();

double time_cutoff = settings::time_cutoff[type().transport_index()];

double distance_cutoff =

(time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;

// Select smaller of the three distances

double distance =

std::min({boundary().distance(), collision_distance(), distance_cutoff});

// Advance particle in space and time

this->move_distance(distance);

double dt = distance / speed;

this->time() += dt;

this->lifetime() += dt;

// Score timed track-length tallies

if (!model::active_timed_tracklength_tallies.empty()) {

score_timed_tracklength_tally(*this, distance);

}

// Score track-length tallies

if (!model::active_tracklength_tallies.empty()) {

score_tracklength_tally(*this, distance);

}

// Score track-length estimate of k-eff

if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {

keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;

}

// Score flux derivative accumulators for differential tallies.

if (!model::active_tallies.empty()) {

score_track_derivative(*this, distance);

}

// Set particle weight to zero if it hit the time boundary

if (distance == distance_cutoff) {

wgt() = 0.0;

}

// Clear surface component if distance is long enough

if (distance > TINY_BIT)

surface() = SURFACE_NONE;

}

void Particle::event_cross_surface()

{

// Saving previous cell data

for (int j = 0; j < n_coord(); ++j) {

cell_last(j) = coord(j).cell();

}

n_coord_last() = n_coord();

// Set surface that particle is on and adjust coordinate levels

surface() = boundary().surface();

n_coord() = boundary().coord_level();

if (boundary().lattice_translation()[0] != 0 ||

boundary().lattice_translation()[1] != 0 ||

boundary().lattice_translation()[2] != 0) {

// Particle crosses lattice boundary

int i_lattice = coord(boundary().coord_level() - 1).lattice();

bool verbose = settings::verbosity >= 10 || trace();

cross_lattice(*this, boundary(), verbose);

event() = TallyEvent::LATTICE;

// Score cell to cell partial currents

if (!model::active_surface_tallies.empty()) {

auto& lat {*model::lattices[i_lattice]};

bool is_valid;

Direction normal =

lat.get_normal(boundary().lattice_translation(), is_valid);

if (is_valid) {

normal /= normal.norm();

score_surface_tally(*this, model::active_surface_tallies, normal);

}

}

} else {

const auto& surf {*model::surfaces[surface_index()].get()};

// Particle crosses surface

// If BC, add particle to surface source before crossing surface

if (surf.surf_source_ && surf.bc_) {

add_surf_source_to_bank(*this, surf);

}

this->cross_surface(surf);

// If no BC, add particle to surface source after crossing surface

if (surf.surf_source_ && !surf.bc_) {

add_surf_source_to_bank(*this, surf);

}

if (settings::weight_window_checkpoint_surface) {

apply_weight_windows(*this);

}

event() = TallyEvent::SURFACE;

// Score cell to cell partial currents

if (!model::active_surface_tallies.empty()) {

Direction normal = surf.normal(r());

normal /= normal.norm();

score_surface_tally(*this, model::active_surface_tallies, normal);

}

}

}

void Particle::event_collide()

{

// Score collision estimate of keff

if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {

keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;

}

// Score surface current tallies -- this has to be done before the collision

// since the direction of the particle will change and we need to use the

// pre-collision direction to figure out what mesh surfaces were crossed

if (!model::active_meshsurf_tallies.empty())

score_meshsurface_tally(*this, model::active_meshsurf_tallies);

// Preserve whether the particle is still associated with a recently crossed

// surface so that a direction change during a near-surface collision can be

// reconciled afterward. The surface marker is no longer needed during the

// collision itself.

const bool near_surface = surface() != SURFACE_NONE;

surface() = SURFACE_NONE;

if (settings::run_CE) {

collision(*this);

} else {

collision_mg(*this);

}

// Collision track feature to recording particle interaction

if (settings::collision_track) {

collision_track_record(*this);

}

// Score collision estimator tallies -- this is done after a collision

// has occurred rather than before because we need information on the

// outgoing energy for any tallies with an outgoing energy filter

if (!model::active_collision_tallies.empty())

score_collision_tally(*this);

if (!model::active_analog_tallies.empty()) {

if (settings::run_CE) {

score_analog_tally_ce(*this);

} else {

score_analog_tally_mg(*this);

}

}

if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {

pht_collision_energy();

}

// Reset banked weight during collision

n_bank() = 0;

bank_second_E() = 0.0;

wgt_bank() = 0.0;

// Clear number of secondaries in this collision. This is

// distinct from the number of created neutrons n_bank() above!

n_secondaries() = 0;

zero_delayed_bank();

// Reset fission logical

fission() = false;

// Save coordinates for tallying purposes

r_last_current() = r();

// Set last material to none since cross sections will need to be

// re-evaluated

material_last() = C_NONE;

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

const auto& u {coord(j).u()};

coord(j + 1).u() = u.rotate(m);

} else {

// Otherwise, copy this level's direction

coord(j + 1).u() = coord(j).u();

}

}

// Score flux derivative accumulators for differential tallies.

if (!model::active_tallies.empty())

score_collision_derivative(*this);

#ifdef OPENMC_DAGMC_ENABLED

history().reset();

#endif

if (near_surface && alive())

reconcile_cell_after_collision(*this);

}

void Particle::event_revive_from_secondary(const SourceSite& site)

{

// Write final position for the previous track (skip if this is a freshly

// constructed particle with no prior track, e.g., Phase 2 of shared

// secondary transport)

if (write_track() && n_event() > 0) {

write_particle_track(*this);

}

from_source(&site);

n_event() = 0;

if (!settings::use_shared_secondary_bank) {

n_tracks()++;

}

bank_second_E() = 0.0;

// Subtract secondary particle energy from interim pulse-height results.

// In shared secondary mode, this subtraction was already done on the parent

// particle during create_secondary(), so skip it here.

if (!settings::use_shared_secondary_bank &&

!model::active_pulse_height_tallies.empty() && this->type().is_photon()) {

// Since the birth cell of the particle has not been set we

// have to determine it before the energy of the secondary particle can be

// removed from the pulse-height of this cell.

if (lowest_coord().cell() == C_NONE) {

bool verbose = settings::verbosity >= 10 || trace();

if (!exhaustive_find_cell(*this, verbose)) {

mark_as_lost("Could not find the cell containing particle " +

std::to_string(id()));

return;

}

// Set birth cell attribute

if (cell_born() == C_NONE)

cell_born() = lowest_coord().cell();

// Initialize last cells from current cell

for (int j = 0; j < n_coord(); ++j) {

cell_last(j) = coord(j).cell();

}

n_coord_last() = n_coord();

}

pht_secondary_particles();

}

// Enter new particle in particle track file

if (write_track())

add_particle_track(*this);

}

void Particle::event_check_limit_and_revive()

{

// If particle has too many events, display warning and kill it

n_event()++;

if (n_event() == settings::max_particle_events) {

warning("Particle " + std::to_string(id()) +

" underwent maximum number of events.");

wgt() = 0.0;

}

// In non-shared-secondary mode, revive from local secondary bank

if (!alive() && !settings::use_shared_secondary_bank &&

!local_secondary_bank().empty()) {

SourceSite& site = local_secondary_bank().back();

event_revive_from_secondary(site);

local_secondary_bank().pop_back();

}

}

void Particle::event_death()

{

#ifdef OPENMC_DAGMC_ENABLED

history().reset();

#endif

// Finish particle track output.

if (write_track()) {

write_particle_track(*this);

finalize_particle_track(*this);

}

// Contribute tally reduction variables to global accumulator

const auto k_absorption = keff_tally_absorption();

const auto k_collision = keff_tally_collision();

const auto k_tracklength = keff_tally_tracklength();

const auto leakage = keff_tally_leakage();

if (settings::run_mode == RunMode::EIGENVALUE) {

if (k_absorption != 0.0) {

#pragma omp atomic

global_tally_absorption += k_absorption;

}

if (k_collision != 0.0) {

#pragma omp atomic

global_tally_collision += k_collision;

}

if (k_tracklength != 0.0) {

#pragma omp atomic

global_tally_tracklength += k_tracklength;

}

}

if (leakage != 0.0) {

#pragma omp atomic

global_tally_leakage += leakage;

}

// Reset particle tallies once accumulated

keff_tally_absorption() = 0.0;

keff_tally_collision() = 0.0;

keff_tally_tracklength() = 0.0;

keff_tally_leakage() = 0.0;

if (!model::active_pulse_height_tallies.empty()) {

score_pulse_height_tally(*this, model::active_pulse_height_tallies);

}

// Accumulate track count for this particle history

if (!settings::use_shared_secondary_bank) {

#pragma omp atomic

simulation::simulation_tracks_completed += n_tracks();

}

// Record the number of progeny created by this particle.

// This data will be used to efficiently sort the fission bank.

if (settings::run_mode == RunMode::EIGENVALUE ||

settings::use_shared_secondary_bank) {

simulation::progeny_per_particle[current_work()] = n_progeny();

}

}

void Particle::pht_collision_energy()

{

// Adds the energy particles lose in a collision to the pulse-height

// determine index of cell in pulse_height_cells

auto it = std::find(model::pulse_height_cells.begin(),

model::pulse_height_cells.end(), lowest_coord().cell());

if (it != model::pulse_height_cells.end()) {

int index = std::distance(model::pulse_height_cells.begin(), it);

pht_storage()[index] += E_last() - E();

// If the energy of the particle is below the cutoff, it will not be sampled

// so its energy is added to the pulse-height in the cell

int photon = ParticleType::photon().transport_index();

if (E() < settings::energy_cutoff[photon]) {

pht_storage()[index] += E();

}

}

}

void Particle::pht_secondary_particles()

{

// Removes the energy of secondary produced particles from the pulse-height

// determine index of cell in pulse_height_cells

auto it = std::find(model::pulse_height_cells.begin(),

model::pulse_height_cells.end(), cell_born());

if (it != model::pulse_height_cells.end()) {

int index = std::distance(model::pulse_height_cells.begin(), it);

pht_storage()[index] -= E();

}

}

void Particle::cross_surface(const Surface& surf)

{

if (settings::verbosity >= 10 || trace()) {

write_message(1, " Crossing surface {}", surf.id_);

}

// Handle any applicable boundary conditions.

if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&

settings::run_mode != RunMode::VOLUME) {

surf.bc_->handle_particle(*this, surf);

return;

}

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

// SEARCH NEIGHBOR LISTS FOR NEXT CELL

#ifdef OPENMC_DAGMC_ENABLED

// in DAGMC, we know what the next cell should be

if (surf.geom_type() == GeometryType::DAG) {

int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),

lowest_coord().universe()) -

1;

// save material, temperature, and density multiplier

material_last() = material();

sqrtkT_last() = sqrtkT();

density_mult_last() = density_mult();

// set new cell value

lowest_coord().cell() = i_cell;

auto& cell = model::cells[i_cell];

cell_instance() = 0;

if (cell->distribcell_index_ >= 0)

cell_instance() = cell_instance_at_level(*this, n_coord() - 1);

material() = cell->material(cell_instance());

sqrtkT() = cell->sqrtkT(cell_instance());

density_mult() = cell->density_mult(cell_instance());

return;

}

#endif

bool verbose = settings::verbosity >= 10 || trace();

if (neighbor_list_find_cell(*this, verbose)) {

return;

}

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

// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS

// Remove lower coordinate levels

n_coord() = 1;

bool found = exhaustive_find_cell(*this, verbose);

if (settings::run_mode != RunMode::PLOTTING && (!found)) {

// If a cell is still not found, there are two possible causes: 1) there is

// a void in the model, and 2) the particle hit a surface at a tangent. If

// the particle is really traveling tangent to a surface, if we move it

// forward a tiny bit it should fix the problem.

surface() = SURFACE_NONE;

n_coord() = 1;

r() += TINY_BIT * u();

// Couldn't find next cell anywhere! This probably means there is an actual

// undefined region in the geometry.

if (!exhaustive_find_cell(*this, verbose)) {

mark_as_lost("After particle " + std::to_string(id()) +

" crossed surface " + std::to_string(surf.id_) +

" it could not be located in any cell and it did not leak.");

return;

}

}

}

void Particle::cross_vacuum_bc(const Surface& surf)

{

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

// still processed

if (!model::active_meshsurf_tallies.empty()) {

// TODO: Find a better solution to score surface currents than

// physically moving the particle forward slightly

r() += TINY_BIT * u();

score_meshsurface_tally(*this, model::active_meshsurf_tallies);

}

// Score to global leakage tally

keff_tally_leakage() += wgt();

// Kill the particle

wgt() = 0.0;

// Display message

if (settings::verbosity >= 10 || trace()) {

write_message(1, " Leaked out of surface {}", surf.id_);

}

}

void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)

{

// Do not handle reflective boundary conditions on lower universes

if (n_coord() != 1) {

mark_as_lost("Cannot reflect particle " + std::to_string(id()) +

" off surface in a lower universe.");

return;

}

// Score surface currents since reflection causes the direction of the

// particle to change. For surface filters, we need to score the tallies

// twice, once before the particle's surface attribute has changed and

// once after. For mesh surface filters, we need to artificially move

// the particle slightly back in case the surface crossing is coincident

// with a mesh boundary

if (!model::active_surface_tallies.empty()) {

Direction normal = surf.normal(r());

normal /= normal.norm();

score_surface_tally(*this, model::active_surface_tallies, normal);

}

if (!model::active_meshsurf_tallies.empty()) {

Position r {this->r()};

this->r() -= TINY_BIT * u();

score_meshsurface_tally(*this, model::active_meshsurf_tallies);

this->r() = r;

}

// Set the new particle direction

u() = new_u;

// Reassign particle's cell and surface

coord(0).cell() = cell_last(0);

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.

// (unless we're using a dagmc model, which has exactly one universe)

n_coord() = 1;

if (surf.geom_type() != GeometryType::DAG &&

!neighbor_list_find_cell(*this)) {

mark_as_lost("Couldn't find particle after reflecting from surface " +

std::to_string(surf.id_) + ".");

return;

}

// Set previous coordinate going slightly past surface crossing

r_last_current() = r() + TINY_BIT * u();

// Diagnostic message

if (settings::verbosity >= 10 || trace()) {

write_message(1, " Reflected from surface {}", surf.id_);

}

}

void Particle::cross_periodic_bc(

const Surface& surf, Position new_r, Direction new_u, int new_surface)

{

// Do not handle periodic boundary conditions on lower universes

if (n_coord() != 1) {

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;

}

// Score surface currents since reflection causes the direction of the

// 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 {this->r()};

this->r() -= TINY_BIT * u();

score_meshsurface_tally(*this, model::active_meshsurf_tallies);

this->r() = r;

}

// Adjust the particle's location and direction.

r() = new_r;

u() = new_u;

// Reassign particle's surface

surface() = new_surface;

// Figure out what cell particle is in now

n_coord() = 1;

if (!neighbor_list_find_cell(*this)) {

mark_as_lost("Couldn't find particle after hitting periodic "

"boundary on surface " +

std::to_string(surf.id_) + ".");

return;

}

// Set previous coordinate going slightly past surface crossing

r_last_current() = r() + TINY_BIT * u();

// Diagnostic message

if (settings::verbosity >= 10 || trace()) {

write_message(1, " Hit periodic boundary on surface {}", surf.id_);

}

}

void Particle::mark_as_lost(const char* message)

{

// Print warning and write lost particle file

warning(message);

if (settings::max_write_lost_particles < 0 ||

simulation::n_lost_particles < settings::max_write_lost_particles) {

write_restart();

}

// Increment number of lost particles

wgt() = 0.0;

#pragma omp atomic

simulation::n_lost_particles += 1;

// Count the total number of simulated particles (on this processor)

auto n = simulation::current_batch * settings::gen_per_batch *

simulation::work_per_rank;

// Abort the simulation if the maximum number of lost particles has been

// reached

if (simulation::n_lost_particles >= settings::max_lost_particles &&

simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {

fatal_error("Maximum number of lost particles has been reached.");

}

}

void Particle::write_restart() const

{

// Dont write another restart file if in particle restart mode

if (settings::run_mode == RunMode::PARTICLE)

return;

// Set up file name

auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,

simulation::current_batch, id());

#pragma omp critical(WriteParticleRestart)

{

// Create file

hid_t file_id = file_open(filename, 'w');

// Write filetype and version info

write_attribute(file_id, "filetype", "particle restart");

write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);

write_attribute(file_id, "openmc_version", VERSION);

#ifdef GIT_SHA1

write_attr_string(file_id, "git_sha1", GIT_SHA1);

#endif

// Write data to file

write_dataset(file_id, "current_batch", simulation::current_batch);

write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);

write_dataset(file_id, "current_generation", simulation::current_gen);

write_dataset(file_id, "n_particles", settings::n_particles);

switch (settings::run_mode) {

case RunMode::FIXED_SOURCE:

write_dataset(file_id, "run_mode", "fixed source");

break;

case RunMode::EIGENVALUE:

write_dataset(file_id, "run_mode", "eigenvalue");

break;

case RunMode::PARTICLE:

write_dataset(file_id, "run_mode", "particle restart");

break;

default:

break;

}

write_dataset(file_id, "id", id());

write_dataset(file_id, "type", type().pdg_number());

// Get source site data for the particle that got lost

int64_t i = current_work();

SourceSite site;

if (settings::run_mode == RunMode::EIGENVALUE) {

site = simulation::source_bank[i];

} else if (settings::run_mode == RunMode::FIXED_SOURCE &&

settings::use_shared_secondary_bank &&

i < simulation::shared_secondary_bank_read.size()) {

site = simulation::shared_secondary_bank_read[i];

} else if (settings::run_mode == RunMode::FIXED_SOURCE) {

// Re-sample using the same seed used to generate the source particle.

// current_work() is 0-indexed, compute_particle_id expects 1-indexed.

int64_t id = compute_transport_seed(compute_particle_id(i + 1));

uint64_t seed = init_seed(id, STREAM_SOURCE);

site = sample_external_source(&seed);

}

write_dataset(file_id, "weight", site.wgt);

write_dataset(file_id, "energy", site.E);

write_dataset(file_id, "xyz", site.r);

write_dataset(file_id, "uvw", site.u);

write_dataset(file_id, "time", site.time);

// Close file

file_close(file_id);

} // #pragma omp critical

}

void Particle::update_neutron_xs(

int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)

{

// Get microscopic cross section cache

auto& micro = this->neutron_xs(i_nuclide);

// If the cache doesn't match, recalculate micro xs

if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||

i_sab != micro.index_sab || sab_frac != micro.sab_frac ||

ncrystal_xs != micro.ncrystal_xs) {

data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);

// If NCrystal is being used, update micro cross section cache

micro.ncrystal_xs = ncrystal_xs;

if (ncrystal_xs >= 0.0) {

data::nuclides[i_nuclide]->calculate_elastic_xs(*this);

ncrystal_update_micro(ncrystal_xs, micro);

}

}

}

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

// Non-method functions

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

void add_surf_source_to_bank(Particle& p, const Surface& surf)

{

if (simulation::current_batch <= settings::n_inactive ||

simulation::surf_source_bank.full()) {

return;

}

// If a cell/cellfrom/cellto parameter is defined

if (settings::ssw_cell_id != C_NONE) {

// Retrieve cell index and storage type

int cell_idx = model::cell_map[settings::ssw_cell_id];