mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Convert transport() and cross_surface() to C++
This commit is contained in:
Paul Romano
parent
5af7510c40
commit
f8925b25be
21 changed files with 566 additions and 700 deletions
|
|
@ -337,7 +337,6 @@ add_library(libopenmc SHARED
|
|||
src/stl_vector.F90
|
||||
src/string.F90
|
||||
src/surface_header.F90
|
||||
src/tracking.F90
|
||||
src/track_output.F90
|
||||
src/vector_header.F90
|
||||
src/xml_interface.F90
|
||||
|
|
|
|||
|
|
@ -195,5 +195,13 @@ public:
|
|||
};
|
||||
#endif
|
||||
|
||||
//==============================================================================
|
||||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
#ifdef DAGMC
|
||||
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed);
|
||||
#endif
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_CELL_H
|
||||
|
|
|
|||
|
|
@ -49,7 +49,7 @@ void warning(const std::stringstream& message)
|
|||
warning(message.str());
|
||||
}
|
||||
|
||||
void write_message(const std::string& message, int level);
|
||||
void write_message(const std::string& message, int level=0);
|
||||
|
||||
inline
|
||||
void write_message(const std::stringstream& message, int level)
|
||||
|
|
|
|||
|
|
@ -148,6 +148,12 @@ extern "C" {
|
|||
//! \param src Source site data
|
||||
void from_source(const Bank* src);
|
||||
|
||||
//! Transport the particle
|
||||
void transport();
|
||||
|
||||
//! Cross a surface
|
||||
void cross_surface();
|
||||
|
||||
//! mark a particle as lost and create a particle restart file
|
||||
//! \param message A warning message to display
|
||||
void mark_as_lost(const char* message);
|
||||
|
|
|
|||
|
|
@ -14,7 +14,7 @@ namespace openmc {
|
|||
//==============================================================================
|
||||
|
||||
//! Sample a nuclide and reaction and then calls the appropriate routine
|
||||
extern "C" void collision(Particle* p);
|
||||
void collision(Particle* p);
|
||||
|
||||
//! Samples an incident neutron reaction
|
||||
void sample_neutron_reaction(Particle* p);
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@ namespace openmc {
|
|||
|
||||
//! \brief samples particle behavior after a collision event.
|
||||
//! \param p Particle to operate on
|
||||
extern "C" void
|
||||
void
|
||||
collision_mg(Particle* p);
|
||||
|
||||
//! \brief samples a reaction type.
|
||||
|
|
|
|||
|
|
@ -35,6 +35,13 @@ void
|
|||
apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
|
||||
double atom_density, int score_bin, double& score);
|
||||
|
||||
void score_collision_derivative(const Particle* p);
|
||||
|
||||
void score_track_derivative(const Particle* p, double distance);
|
||||
|
||||
//! Set the flux derivatives on differential tallies to zero.
|
||||
void zero_flux_derivs();
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -50,6 +50,16 @@ private:
|
|||
// Non-member functions
|
||||
//==============================================================================
|
||||
|
||||
void score_collision_tally(const Particle* p);
|
||||
|
||||
void score_analog_tally_ce(const Particle* p);
|
||||
|
||||
void score_analog_tally_mg(const Particle* p);
|
||||
|
||||
void score_tracklength_tally(const Particle* p, double distance);
|
||||
|
||||
void score_surface_tally(const Particle* p, const std::vector<int>& tallies);
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_TALLIES_TALLY_SCORING_H
|
||||
|
|
|
|||
|
|
@ -805,12 +805,12 @@ extern "C" {
|
|||
}
|
||||
|
||||
#ifdef DAGMC
|
||||
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed )
|
||||
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed)
|
||||
{
|
||||
moab::EntityHandle surf =
|
||||
surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_);
|
||||
surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_);
|
||||
moab::EntityHandle vol =
|
||||
cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_);
|
||||
cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_);
|
||||
|
||||
moab::EntityHandle new_vol;
|
||||
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
|
||||
|
|
|
|||
|
|
@ -50,34 +50,10 @@ module geometry
|
|||
integer(C_INT), intent(out) :: lattice_translation(3)
|
||||
integer(C_INT), intent(out) :: next_level
|
||||
end subroutine distance_to_boundary
|
||||
|
||||
#ifdef DAGMC
|
||||
|
||||
function next_cell_c(current_cell, surface_crossed) &
|
||||
bind(C, name="next_cell") result(new_cell)
|
||||
import C_PTR, C_INT32_T
|
||||
type(C_PTR), intent(in), value :: current_cell
|
||||
type(C_PTR), intent(in), value :: surface_crossed
|
||||
integer(C_INT32_T) :: new_cell
|
||||
end function next_cell_c
|
||||
|
||||
#endif
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
#ifdef DAGMC
|
||||
|
||||
function next_cell(c, s) result(new_cell)
|
||||
type(Cell), intent(in) :: c
|
||||
type(Surface), intent(in) :: s
|
||||
integer :: new_cell
|
||||
new_cell = next_cell_c(c%ptr, s%ptr)
|
||||
end function next_cell
|
||||
|
||||
#endif
|
||||
|
||||
!===============================================================================
|
||||
! FIND_CELL determines what cell a source particle is in within a particular
|
||||
! universe. If the base universe is passed, the particle should be found as long
|
||||
|
|
|
|||
|
|
@ -1135,11 +1135,6 @@ extern "C" {
|
|||
return model::materials[i_mat - 1]->density_gpcc_;
|
||||
}
|
||||
|
||||
void material_calculate_xs(const Particle* p)
|
||||
{
|
||||
model::materials[p->material - 1]->calculate_xs(*p);
|
||||
}
|
||||
|
||||
void free_memory_material()
|
||||
{
|
||||
for (Material *mat : model::materials) {delete mat;}
|
||||
|
|
|
|||
|
|
@ -2,12 +2,9 @@ module material_header
|
|||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use particle_header, only: Particle
|
||||
|
||||
implicit none
|
||||
|
||||
private
|
||||
public :: material_calculate_xs
|
||||
public :: material_id
|
||||
public :: material_nuclide
|
||||
public :: material_nuclide_size
|
||||
|
|
@ -21,11 +18,6 @@ module material_header
|
|||
integer(C_INT32_T) :: id
|
||||
end function
|
||||
|
||||
subroutine material_calculate_xs(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(in) :: p
|
||||
end subroutine
|
||||
|
||||
function material_nuclide(i_mat, idx) bind(C) result(nuc)
|
||||
import C_INT32_T, C_INT
|
||||
integer(C_INT32_T), value :: i_mat
|
||||
|
|
|
|||
509
src/particle.cpp
509
src/particle.cpp
|
|
@ -1,16 +1,29 @@
|
|||
#include "openmc/particle.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <algorithm> // copy, min
|
||||
#include <cmath> // log, abs, copysign
|
||||
#include <sstream>
|
||||
|
||||
#include "openmc/bank.h"
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/cell.h"
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/error.h"
|
||||
#include "openmc/geometry.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/material.h"
|
||||
#include "openmc/message_passing.h"
|
||||
#include "openmc/mgxs_interface.h"
|
||||
#include "openmc/nuclide.h"
|
||||
#include "openmc/physics.h"
|
||||
#include "openmc/physics_mg.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/surface.h"
|
||||
#include "openmc/simulation.h"
|
||||
#include "openmc/tallies/derivative.h"
|
||||
#include "openmc/tallies/tally.h"
|
||||
#include "openmc/tallies/tally_scoring.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
|
|
@ -119,6 +132,500 @@ Particle::from_source(const Bank* src)
|
|||
last_E = E;
|
||||
}
|
||||
|
||||
extern "C" void initialize_particle_track();
|
||||
extern "C" void write_particle_track(const Particle*);
|
||||
extern "C" void add_particle_track();
|
||||
extern "C" void finalize_particle_track(const Particle*);
|
||||
|
||||
void
|
||||
Particle::transport()
|
||||
{
|
||||
// Display message if high verbosity or trace is on
|
||||
if (settings::verbosity >= 9 || simulation::trace) {
|
||||
write_message("Simulating Particle " + std::to_string(id));
|
||||
}
|
||||
|
||||
// Initialize number of events to zero
|
||||
int n_event = 0;
|
||||
|
||||
// Add paricle's starting weight to count for normalizing tallies later
|
||||
#pragma omp atomic
|
||||
total_weight += wgt;
|
||||
|
||||
// Force calculation of cross-sections by setting last energy to zero
|
||||
if (settings::run_CE) {
|
||||
for (int i = 0; i < data::nuclides.size(); ++i) {
|
||||
simulation::micro_xs[i].last_E = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Prepare to write out particle track.
|
||||
if (write_track) initialize_particle_track();
|
||||
|
||||
// Every particle starts with no accumulated flux derivative.
|
||||
if (!model::active_tallies.empty()) zero_flux_derivs();
|
||||
|
||||
while (true) {
|
||||
// Set the random number stream
|
||||
if (type == static_cast<int>(ParticleType::neutron)) {
|
||||
prn_set_stream(STREAM_TRACKING);
|
||||
} else {
|
||||
prn_set_stream(STREAM_PHOTON);
|
||||
}
|
||||
|
||||
// Store pre-collision particle properties
|
||||
last_wgt = wgt;
|
||||
last_E = E;
|
||||
std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw);
|
||||
std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz);
|
||||
|
||||
// If the cell hasn't been determined based on the particle's location,
|
||||
// initiate a search for the current cell. This generally happens at the
|
||||
// beginning of the history and again for any secondary particles
|
||||
if (coord[n_coord - 1].cell == C_NONE) {
|
||||
if (!find_cell(this, false)) {
|
||||
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 = coord[n_coord - 1].cell;
|
||||
}
|
||||
|
||||
// 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 != last_material || sqrtkT != last_sqrtkT) {
|
||||
// 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 - 1]->calculate_xs(*this);
|
||||
}
|
||||
} else {
|
||||
// Get the MG data
|
||||
calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw,
|
||||
simulation::material_xs.total, simulation::material_xs.absorption,
|
||||
simulation::material_xs.nu_fission);
|
||||
|
||||
// Finally, update the particle group while we have already checked
|
||||
// for if multi-group
|
||||
last_g = g;
|
||||
}
|
||||
} else {
|
||||
simulation::material_xs.total = 0.0;
|
||||
simulation::material_xs.absorption = 0.0;
|
||||
simulation::material_xs.fission = 0.0;
|
||||
simulation::material_xs.nu_fission = 0.0;
|
||||
}
|
||||
|
||||
// Find the distance to the nearest boundary
|
||||
double d_boundary;
|
||||
int surface_crossed;
|
||||
int lattice_translation[3];
|
||||
int next_level;
|
||||
distance_to_boundary(this, &d_boundary, &surface_crossed,
|
||||
lattice_translation, &next_level);
|
||||
|
||||
// Sample a distance to collision
|
||||
double d_collision;
|
||||
if (type == static_cast<int>(ParticleType::electron) ||
|
||||
type == static_cast<int>(ParticleType::positron)) {
|
||||
d_collision = 0.0;
|
||||
} else if (simulation::material_xs.total == 0.0) {
|
||||
d_collision = INFINITY;
|
||||
} else {
|
||||
d_collision = -std::log(prn()) / simulation::material_xs.total;
|
||||
}
|
||||
|
||||
// Select smaller of the two distances
|
||||
double distance = std::min(d_boundary, d_collision);
|
||||
|
||||
// Advance particle
|
||||
for (int j = 0; j < n_coord; ++j) {
|
||||
// TODO: use Position
|
||||
coord[j].xyz[0] += distance * coord[j].uvw[0];
|
||||
coord[j].xyz[1] += distance * coord[j].uvw[1];
|
||||
coord[j].xyz[2] += distance * coord[j].uvw[2];
|
||||
}
|
||||
|
||||
// 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 == RUN_MODE_EIGENVALUE &&
|
||||
type == static_cast<int>(ParticleType::neutron)) {
|
||||
global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission;
|
||||
}
|
||||
|
||||
// Score flux derivative accumulators for differential tallies.
|
||||
if (!model::active_tallies.empty()) {
|
||||
score_track_derivative(this, distance);
|
||||
}
|
||||
|
||||
if (d_collision > d_boundary) {
|
||||
// ====================================================================
|
||||
// PARTICLE CROSSES SURFACE
|
||||
|
||||
if (next_level > 0) n_coord = next_level;
|
||||
|
||||
// Saving previous cell data
|
||||
for (int j = 0; j < n_coord; ++j) {
|
||||
last_cell[j] = coord[j].cell;
|
||||
}
|
||||
last_n_coord = n_coord;
|
||||
|
||||
if (lattice_translation[0] != 0 || lattice_translation[1] != 0 ||
|
||||
lattice_translation[2] != 0) {
|
||||
// Particle crosses lattice boundary
|
||||
surface = ERROR_INT;
|
||||
cross_lattice(this, lattice_translation);
|
||||
event = EVENT_LATTICE;
|
||||
} else {
|
||||
// Particle crosses surface
|
||||
surface = surface_crossed;
|
||||
this->cross_surface();
|
||||
event = EVENT_SURFACE;
|
||||
}
|
||||
// Score cell to cell partial currents
|
||||
if (!model::active_surface_tallies.empty()) {
|
||||
score_surface_tally(this, model::active_surface_tallies);
|
||||
}
|
||||
} else {
|
||||
// ====================================================================
|
||||
// PARTICLE HAS COLLISION
|
||||
|
||||
// Score collision estimate of keff
|
||||
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
|
||||
type == static_cast<int>(ParticleType::neutron)) {
|
||||
global_tally_collision += wgt * simulation::material_xs.nu_fission
|
||||
/ simulation::material_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_surface_tally(this, model::active_meshsurf_tallies);
|
||||
|
||||
// Clear surface component
|
||||
surface = ERROR_INT;
|
||||
|
||||
if (settings::run_CE) {
|
||||
collision(this);
|
||||
} else {
|
||||
collision_mg(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);
|
||||
}
|
||||
}
|
||||
|
||||
// Reset banked weight during collision
|
||||
n_bank = 0;
|
||||
wgt_bank = 0.0;
|
||||
for (int& v : n_delayed_bank) v = 0;
|
||||
|
||||
// Reset fission logical
|
||||
fission = false;
|
||||
|
||||
// Save coordinates for tallying purposes
|
||||
std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current);
|
||||
|
||||
// Set last material to none since cross sections will need to be
|
||||
// re-evaluated
|
||||
last_material = F90_NONE;
|
||||
|
||||
// Set all uvws to base level -- right now, after a collision, only the
|
||||
// base level uvws are changed
|
||||
for (int j = 0; j < n_coord - 1; ++j) {
|
||||
if (coord[j + 1].rotated) {
|
||||
// If next level is rotated, apply rotation matrix
|
||||
// FIXME
|
||||
// coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % &
|
||||
// rotation_matrix, coord[j].uvw)
|
||||
} else {
|
||||
// Otherwise, copy this level's direction
|
||||
std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw);
|
||||
}
|
||||
}
|
||||
|
||||
// Score flux derivative accumulators for differential tallies.
|
||||
if (!model::active_tallies.empty()) score_collision_derivative(this);
|
||||
}
|
||||
|
||||
// If particle has too many events, display warning and kill it
|
||||
++n_event;
|
||||
if (n_event == MAX_EVENTS) {
|
||||
if (mpi::master) warning("Particle " + std::to_string(id)
|
||||
+ " underwent maximum number of events.");
|
||||
alive = false;
|
||||
}
|
||||
|
||||
// Check for secondary particles if this particle is dead
|
||||
if (!alive) {
|
||||
// If no secondary particles, break out of event loop
|
||||
if (n_secondary == 0) break;
|
||||
|
||||
this->from_source(&secondary_bank[n_secondary - 1]);
|
||||
--n_secondary;
|
||||
n_event = 0;
|
||||
|
||||
// Enter new particle in particle track file
|
||||
if (write_track) add_particle_track();
|
||||
}
|
||||
}
|
||||
|
||||
// Finish particle track output.
|
||||
if (write_track) {
|
||||
write_particle_track(this);
|
||||
finalize_particle_track(this);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
Particle::cross_surface()
|
||||
{
|
||||
int i_surface = std::abs(surface);
|
||||
// TODO: off-by-one
|
||||
const auto& surf {model::surfaces[i_surface - 1]};
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
write_message(" Crossing surface " + std::to_string(surf->id_));
|
||||
}
|
||||
|
||||
if (surf->bc_ == BC_VACUUM && (settings::run_mode != RUN_MODE_PLOTTING)) {
|
||||
// =======================================================================
|
||||
// PARTICLE LEAKS OUT OF PROBLEM
|
||||
|
||||
// Kill particle
|
||||
alive = false;
|
||||
|
||||
// Score any surface current tallies -- note that the particle is moved
|
||||
// forward slightly so that if the mesh boundary is on the surface, it is
|
||||
// still processed
|
||||
|
||||
if (!model::active_meshsurf_tallies.empty()) {
|
||||
// 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];
|
||||
score_surface_tally(this, model::active_meshsurf_tallies);
|
||||
}
|
||||
|
||||
// Score to global leakage tally
|
||||
global_tally_leakage += wgt;
|
||||
|
||||
// Display message
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
write_message(" Leaked out of surface " + std::to_string(surf->id_));
|
||||
}
|
||||
return;
|
||||
|
||||
} else if (surf->bc_ == BC_REFLECT && (settings::run_mode != RUN_MODE_PLOTTING)) {
|
||||
// =======================================================================
|
||||
// PARTICLE REFLECTS FROM SURFACE
|
||||
|
||||
// Do not handle reflective boundary conditions on lower universes
|
||||
if (n_coord != 1) {
|
||||
this->mark_as_lost("Cannot reflect particle " + std::to_string(id) +
|
||||
" 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()) {
|
||||
score_surface_tally(this, model::active_surface_tallies);
|
||||
}
|
||||
|
||||
|
||||
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];
|
||||
score_surface_tally(this, model::active_meshsurf_tallies);
|
||||
std::copy(&r.x, &r.x + 3, coord[0].xyz);
|
||||
}
|
||||
|
||||
// Reflect particle off surface
|
||||
Direction u = surf->reflect(coord[0].xyz, coord[0].uvw);
|
||||
|
||||
// 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;
|
||||
|
||||
// Reassign particle's cell and surface
|
||||
coord[0].cell = last_cell[last_n_coord - 1];
|
||||
surface = -surface;
|
||||
|
||||
// If a reflective surface is coincident with a lattice or universe
|
||||
// boundary, it is necessary to redetermine the particle's coordinates in
|
||||
// the lower universes.
|
||||
|
||||
n_coord = 1;
|
||||
if (!find_cell(this, true)) {
|
||||
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
|
||||
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];
|
||||
|
||||
// Diagnostic message
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
write_message(" Reflected from surface " + std::to_string(surf->id_));
|
||||
}
|
||||
return;
|
||||
|
||||
} else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) {
|
||||
// =======================================================================
|
||||
// PERIODIC BOUNDARY
|
||||
|
||||
// Do not handle periodic boundary conditions on lower universes
|
||||
if (n_coord != 1) {
|
||||
this->mark_as_lost("Cannot transfer particle " + std::to_string(id) +
|
||||
" 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 {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];
|
||||
score_surface_tally(this, model::active_meshsurf_tallies);
|
||||
std::copy(&r.x, &r.x + 3, coord[0].xyz);
|
||||
}
|
||||
|
||||
// Get a pointer to the partner periodic surface. Offset the index to
|
||||
// correct for C vs. Fortran indexing.
|
||||
auto surf_p = dynamic_cast<PeriodicSurface*>(surf);
|
||||
if (surf_p) {
|
||||
auto other = dynamic_cast<PeriodicSurface*>(
|
||||
model::surfaces[surf_p->i_periodic_]);
|
||||
|
||||
// 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);
|
||||
|
||||
// Reassign particle's surface
|
||||
// TODO: off-by-one
|
||||
surface = rotational ?
|
||||
surf_p->i_periodic_ + 1 :
|
||||
std::copysign(surf_p->i_periodic_ + 1, surface);
|
||||
}
|
||||
|
||||
// Figure out what cell particle is in now
|
||||
n_coord = 1;
|
||||
|
||||
if (!find_cell(this, true)) {
|
||||
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
|
||||
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];
|
||||
|
||||
// Diagnostic message
|
||||
if (settings::verbosity >= 10 || simulation::trace) {
|
||||
write_message(" Hit periodic boundary on surface " +
|
||||
std::to_string(surf->id_));
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
// ==========================================================================
|
||||
// SEARCH NEIGHBOR LISTS FOR NEXT CELL
|
||||
|
||||
#ifdef DAGMC
|
||||
if (settings::dagmc) {
|
||||
int32_t i_cell = next_cell(model::cells[last_cell[0]],
|
||||
model::surfaces[std::abs(surface)]);
|
||||
// save material and temp
|
||||
last_material = material;
|
||||
last_sqrtkT = sqrtKT;
|
||||
// set new cell value
|
||||
coord[0].cell = i_cell
|
||||
cell_instance = 0;
|
||||
material = model::cells[i_cell]->material[0];
|
||||
sqrtKT = model::cells[i_cell]->sqrtKT[0];
|
||||
return
|
||||
}
|
||||
#endif
|
||||
|
||||
if (find_cell(this, true)) return;
|
||||
|
||||
// ==========================================================================
|
||||
// COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
|
||||
|
||||
// Remove lower coordinate levels and assignment of surface
|
||||
surface = ERROR_INT;
|
||||
n_coord = 1;
|
||||
bool found = find_cell(this, false);
|
||||
|
||||
if (settings::run_mode != RUN_MODE_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.
|
||||
|
||||
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];
|
||||
|
||||
// Couldn't find next cell anywhere// This probably means there is an actual
|
||||
// undefined region in the geometry.
|
||||
|
||||
if (!find_cell(this, false)) {
|
||||
this->mark_as_lost("After particle " + std::to_string(id) +
|
||||
" crossed surface " + std::to_string(surf->id_) +
|
||||
" it could not be located in any cell and it did not leak.");
|
||||
return;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
Particle::mark_as_lost(const char* message)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -74,7 +74,8 @@ contains
|
|||
call set_particle_seed(particle_seed)
|
||||
|
||||
! Transport neutron
|
||||
call transport(p)
|
||||
! FIXME
|
||||
!call transport(p)
|
||||
|
||||
! Write output if particle made it
|
||||
call print_particle(p)
|
||||
|
|
|
|||
|
|
@ -29,11 +29,9 @@ namespace openmc {
|
|||
// data/functions from Fortran side
|
||||
extern "C" void accumulate_tallies();
|
||||
extern "C" void allocate_tally_results();
|
||||
extern "C" void init_tally_routines();
|
||||
extern "C" void setup_active_tallies();
|
||||
extern "C" void simulation_init_f();
|
||||
extern "C" void simulation_finalize_f();
|
||||
extern "C" void transport(Particle* p);
|
||||
extern "C" void write_tallies();
|
||||
|
||||
} // namespace openmc
|
||||
|
|
@ -76,9 +74,6 @@ int openmc_simulation_init()
|
|||
simulation::filter_matches.resize(model::tally_filters.size());
|
||||
}
|
||||
|
||||
// Set up tally procedure pointers
|
||||
init_tally_routines();
|
||||
|
||||
// Allocate source bank, and for eigenvalue simulations also allocate the
|
||||
// fission bank
|
||||
allocate_banks();
|
||||
|
|
@ -213,7 +208,7 @@ int openmc_next_batch(int* status)
|
|||
initialize_history(&p, simulation::current_work);
|
||||
|
||||
// transport particle
|
||||
transport(&p);
|
||||
p.transport();
|
||||
}
|
||||
|
||||
// Accumulate time for transport
|
||||
|
|
|
|||
|
|
@ -340,7 +340,7 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other,
|
|||
Position& r, Direction& u) const
|
||||
{
|
||||
Direction other_n = other->normal(r);
|
||||
if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) {
|
||||
if (other_n.x == 1 && other_n.y == 0 && other_n.z == 0) {
|
||||
r.x = x0_;
|
||||
return false;
|
||||
} else {
|
||||
|
|
@ -397,11 +397,11 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
|
|||
write_dataset(group_id, "coefficients", coeffs);
|
||||
}
|
||||
|
||||
bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other,
|
||||
bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other,
|
||||
Position& r, Direction& u) const
|
||||
{
|
||||
Direction other_n = other->normal(r);
|
||||
if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) {
|
||||
if (other_n.x == 0 && other_n.y == 1 && other_n.z == 0) {
|
||||
// The periodic partner is also aligned along y. Just change the y coord.
|
||||
r.y = y0_;
|
||||
return false;
|
||||
|
|
|
|||
|
|
@ -571,7 +571,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
|
|||
|
||||
//! Adjust diff tally flux derivatives for a particle tracking event.
|
||||
|
||||
extern "C" void
|
||||
void
|
||||
score_track_derivative(const Particle* p, double distance)
|
||||
{
|
||||
// A void material cannot be perturbed so it will not affect flux derivatives.
|
||||
|
|
@ -627,8 +627,7 @@ score_track_derivative(const Particle* p, double distance)
|
|||
//! absorption will never be tallied. The paricle will be killed before any
|
||||
//! further tallies are scored.
|
||||
|
||||
extern "C" void
|
||||
score_collision_derivative(const Particle* p)
|
||||
void score_collision_derivative(const Particle* p)
|
||||
{
|
||||
// A void material cannot be perturbed so it will not affect flux derivatives.
|
||||
if (p->material == MATERIAL_VOID) return;
|
||||
|
|
@ -698,9 +697,7 @@ score_collision_derivative(const Particle* p)
|
|||
}
|
||||
}
|
||||
|
||||
//! Set the flux derivatives on differential tallies to zero.
|
||||
extern "C" void
|
||||
zero_flux_derivs()
|
||||
void zero_flux_derivs()
|
||||
{
|
||||
for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -11,87 +11,16 @@ module tally
|
|||
use message_passing
|
||||
use mgxs_interface
|
||||
use nuclide_header
|
||||
use particle_header, only: LocalCoord, Particle
|
||||
use settings
|
||||
use simulation_header
|
||||
use string, only: to_str
|
||||
use tally_derivative_header
|
||||
use tally_filter
|
||||
use tally_header
|
||||
|
||||
implicit none
|
||||
|
||||
procedure(score_analog_tally_), pointer :: score_analog_tally => null()
|
||||
|
||||
abstract interface
|
||||
subroutine score_analog_tally_(p)
|
||||
import Particle
|
||||
type(Particle), intent(in) :: p
|
||||
end subroutine score_analog_tally_
|
||||
end interface
|
||||
|
||||
interface
|
||||
subroutine score_analog_tally_ce(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(in) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine score_analog_tally_mg(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(in) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine score_tracklength_tally(p, distance) bind(C)
|
||||
import Particle, C_DOUBLE
|
||||
type(Particle) :: p
|
||||
real(C_DOUBLE), value :: distance
|
||||
end subroutine
|
||||
|
||||
subroutine score_collision_tally(p) bind(C)
|
||||
import Particle
|
||||
type(Particle) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine score_meshsurface_tally(p) bind(C)
|
||||
import Particle
|
||||
type(Particle) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine score_surface_tally(p) bind(C)
|
||||
import Particle
|
||||
type(Particle) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine score_track_derivative(p, distance) bind(C)
|
||||
import Particle, C_DOUBLE
|
||||
type(Particle) :: p
|
||||
real(C_DOUBLE), value :: distance
|
||||
end subroutine
|
||||
|
||||
subroutine score_collision_derivative(p) bind(C)
|
||||
import Particle
|
||||
type(Particle) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine zero_flux_derivs() bind(C)
|
||||
end subroutine
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! INIT_TALLY_ROUTINES Sets the procedure pointers needed for minimizing code
|
||||
! with the CE and MG modes.
|
||||
!===============================================================================
|
||||
|
||||
subroutine init_tally_routines() bind(C)
|
||||
if (run_CE) then
|
||||
score_analog_tally => score_analog_tally_ce
|
||||
else
|
||||
score_analog_tally => score_analog_tally_mg
|
||||
end if
|
||||
end subroutine init_tally_routines
|
||||
|
||||
!===============================================================================
|
||||
! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history
|
||||
! within the batch to a new random variable
|
||||
|
|
|
|||
|
|
@ -1976,8 +1976,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux,
|
|||
//
|
||||
//! Analog tallies ar etriggered at every collision, not every event.
|
||||
|
||||
extern "C" void
|
||||
score_analog_tally_ce(const Particle* p)
|
||||
void score_analog_tally_ce(const Particle* p)
|
||||
{
|
||||
for (auto i_tally : model::active_analog_tallies) {
|
||||
const Tally& tally {*model::tallies[i_tally]};
|
||||
|
|
@ -2040,8 +2039,7 @@ score_analog_tally_ce(const Particle* p)
|
|||
//
|
||||
//! Analog tallies ar etriggered at every collision, not every event.
|
||||
|
||||
extern "C" void
|
||||
score_analog_tally_mg(const Particle* p)
|
||||
void score_analog_tally_mg(const Particle* p)
|
||||
{
|
||||
for (auto i_tally : model::active_analog_tallies) {
|
||||
const Tally& tally {*model::tallies[i_tally]};
|
||||
|
|
@ -2096,7 +2094,7 @@ score_analog_tally_mg(const Particle* p)
|
|||
//! collision) and thus cannot be done for tallies that require post-collision
|
||||
//! information.
|
||||
|
||||
extern "C" void
|
||||
void
|
||||
score_tracklength_tally(const Particle* p, double distance)
|
||||
{
|
||||
// Determine the tracklength estimate of the flux
|
||||
|
|
@ -2172,8 +2170,7 @@ score_tracklength_tally(const Particle* p, double distance)
|
|||
//! reactions that we didn't sample. It is assumed the material is not void
|
||||
//! since collisions do not occur in voids.
|
||||
|
||||
extern "C" void
|
||||
score_collision_tally(const Particle* p)
|
||||
void score_collision_tally(const Particle* p)
|
||||
{
|
||||
// Determine the collision estimate of the flux
|
||||
double flux;
|
||||
|
|
@ -2243,8 +2240,8 @@ score_collision_tally(const Particle* p)
|
|||
|
||||
//! Score surface or mesh-surface tallies for particle currents.
|
||||
|
||||
static void
|
||||
score_surface_tally_inner(const Particle* p, const std::vector<int>& tallies)
|
||||
void
|
||||
score_surface_tally(const Particle* p, const std::vector<int>& tallies)
|
||||
{
|
||||
// No collision, so no weight change when survival biasing
|
||||
double flux = p->wgt;
|
||||
|
|
@ -2289,20 +2286,4 @@ score_surface_tally_inner(const Particle* p, const std::vector<int>& tallies)
|
|||
match.bins_present_ = false;
|
||||
}
|
||||
|
||||
//! Score mesh-surface tallies for particle currents.
|
||||
|
||||
extern "C" void
|
||||
score_meshsurface_tally(const Particle* p)
|
||||
{
|
||||
score_surface_tally_inner(p, model::active_meshsurf_tallies);
|
||||
}
|
||||
|
||||
//! Score surface tallies for particle currents.
|
||||
|
||||
extern "C" void
|
||||
score_surface_tally(const Particle* p)
|
||||
{
|
||||
score_surface_tally_inner(p, model::active_surface_tallies);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -5,6 +5,8 @@
|
|||
|
||||
module track_output
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants
|
||||
use hdf5_interface
|
||||
use particle_header, only: Particle
|
||||
|
|
@ -34,7 +36,7 @@ contains
|
|||
! information
|
||||
!===============================================================================
|
||||
|
||||
subroutine initialize_particle_track()
|
||||
subroutine initialize_particle_track() bind(C)
|
||||
allocate(tracks(1))
|
||||
end subroutine initialize_particle_track
|
||||
|
||||
|
|
@ -42,7 +44,7 @@ contains
|
|||
! WRITE_PARTICLE_TRACK copies particle position to an array.
|
||||
!===============================================================================
|
||||
|
||||
subroutine write_particle_track(p)
|
||||
subroutine write_particle_track(p) bind(C)
|
||||
type(Particle), intent(in) :: p
|
||||
real(8), allocatable :: new_coords(:, :)
|
||||
|
||||
|
|
@ -71,7 +73,7 @@ contains
|
|||
! secondary particle
|
||||
!===============================================================================
|
||||
|
||||
subroutine add_particle_track()
|
||||
subroutine add_particle_track() bind(C)
|
||||
type(TrackCoordinates), allocatable :: new_tracks(:)
|
||||
|
||||
integer :: i
|
||||
|
|
@ -92,7 +94,7 @@ contains
|
|||
! FINALIZE_PARTICLE_TRACK writes the particle track array to disk.
|
||||
!===============================================================================
|
||||
|
||||
subroutine finalize_particle_track(p)
|
||||
subroutine finalize_particle_track(p) bind(C)
|
||||
type(Particle), intent(in) :: p
|
||||
|
||||
integer :: i
|
||||
|
|
|
|||
539
src/tracking.F90
539
src/tracking.F90
|
|
@ -1,539 +0,0 @@
|
|||
module tracking
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use constants
|
||||
use error, only: warning, write_message
|
||||
use geometry_header, only: cells
|
||||
use geometry, only: find_cell, distance_to_boundary, cross_lattice,&
|
||||
check_cell_overlap
|
||||
#ifdef DAGMC
|
||||
use geometry, only: next_cell
|
||||
#endif
|
||||
|
||||
use material_header, only: material_calculate_xs
|
||||
use message_passing
|
||||
use mgxs_interface
|
||||
use nuclide_header
|
||||
use particle_header
|
||||
use random_lcg, only: prn, prn_set_stream
|
||||
use settings
|
||||
use simulation_header
|
||||
use string, only: to_str
|
||||
use surface_header
|
||||
use tally_header
|
||||
use tally, only: score_analog_tally, score_tracklength_tally, &
|
||||
score_collision_tally, score_surface_tally, &
|
||||
score_meshsurface_tally, &
|
||||
score_track_derivative, zero_flux_derivs, &
|
||||
score_collision_derivative
|
||||
use track_output, only: initialize_particle_track, write_particle_track, &
|
||||
add_particle_track, finalize_particle_track
|
||||
|
||||
implicit none
|
||||
|
||||
interface
|
||||
subroutine collision(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine
|
||||
|
||||
subroutine collision_mg(p) bind(C)
|
||||
import Particle, C_DOUBLE
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine collision_mg
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! TRANSPORT encompasses the main logic for moving a particle through geometry.
|
||||
!===============================================================================
|
||||
|
||||
subroutine transport(p) bind(C)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
integer :: j ! coordinate level
|
||||
integer :: next_level ! next coordinate level to check
|
||||
integer :: surface_crossed ! surface which particle is on
|
||||
integer :: lattice_translation(3) ! in-lattice translation vector
|
||||
integer :: n_event ! number of collisions/crossings
|
||||
real(8) :: d_boundary ! distance to nearest boundary
|
||||
real(8) :: d_collision ! sampled distance to collision
|
||||
real(8) :: distance ! distance particle travels
|
||||
logical :: found_cell ! found cell which particle is in?
|
||||
|
||||
! Display message if high verbosity or trace is on
|
||||
if (verbosity >= 9 .or. trace) then
|
||||
call write_message("Simulating Particle " // trim(to_str(p % id)))
|
||||
end if
|
||||
|
||||
! Initialize number of events to zero
|
||||
n_event = 0
|
||||
|
||||
! Add paricle's starting weight to count for normalizing tallies later
|
||||
!$omp atomic
|
||||
total_weight = total_weight + p % wgt
|
||||
|
||||
! Force calculation of cross-sections by setting last energy to zero
|
||||
if (run_CE) then
|
||||
micro_xs % last_E = ZERO
|
||||
end if
|
||||
|
||||
! Prepare to write out particle track.
|
||||
if (p % write_track) then
|
||||
call initialize_particle_track()
|
||||
endif
|
||||
|
||||
! Every particle starts with no accumulated flux derivative.
|
||||
if (active_tallies_size() > 0) call zero_flux_derivs()
|
||||
|
||||
EVENT_LOOP: do
|
||||
! Set the random number stream
|
||||
if (p % type == NEUTRON) then
|
||||
call prn_set_stream(STREAM_TRACKING)
|
||||
else
|
||||
call prn_set_stream(STREAM_PHOTON)
|
||||
end if
|
||||
|
||||
! Store pre-collision particle properties
|
||||
p % last_wgt = p % wgt
|
||||
p % last_E = p % E
|
||||
p % last_uvw = p % coord(1) % uvw
|
||||
p % last_xyz = p % coord(1) % xyz
|
||||
|
||||
! 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 (p % coord(p % n_coord) % cell == C_NONE) then
|
||||
call find_cell(p, found_cell)
|
||||
if (.not. found_cell) then
|
||||
call particle_mark_as_lost(p, "Could not find the cell containing" &
|
||||
// " particle " // trim(to_str(p %id)))
|
||||
return
|
||||
end if
|
||||
|
||||
! set birth cell attribute
|
||||
if (p % cell_born == C_NONE) p % cell_born = p % coord(p % n_coord) % cell
|
||||
end if
|
||||
|
||||
! Write particle track.
|
||||
if (p % write_track) call write_particle_track(p)
|
||||
|
||||
if (check_overlaps) call check_cell_overlap(p)
|
||||
|
||||
! Calculate microscopic and macroscopic cross sections
|
||||
if (p % material /= MATERIAL_VOID) then
|
||||
if (run_CE) then
|
||||
if (p % material /= p % last_material .or. &
|
||||
p % sqrtkT /= p % last_sqrtkT) then
|
||||
! 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.
|
||||
call material_calculate_xs(p)
|
||||
end if
|
||||
else
|
||||
! Get the MG data
|
||||
call calculate_xs_c(p % material, p % g, p % sqrtkT, &
|
||||
p % coord(p % n_coord) % uvw, material_xs % total, &
|
||||
material_xs % absorption, material_xs % nu_fission)
|
||||
|
||||
|
||||
! Finally, update the particle group while we have already checked
|
||||
! for if multi-group
|
||||
p % last_g = p % g
|
||||
end if
|
||||
else
|
||||
material_xs % total = ZERO
|
||||
material_xs % absorption = ZERO
|
||||
material_xs % fission = ZERO
|
||||
material_xs % nu_fission = ZERO
|
||||
end if
|
||||
|
||||
! Find the distance to the nearest boundary
|
||||
call distance_to_boundary(p, d_boundary, surface_crossed, &
|
||||
lattice_translation, next_level)
|
||||
|
||||
! Sample a distance to collision
|
||||
if (p % type == ELECTRON .or. p % type == POSITRON) then
|
||||
d_collision = ZERO
|
||||
else if (material_xs % total == ZERO) then
|
||||
d_collision = INFINITY
|
||||
else
|
||||
d_collision = -log(prn()) / material_xs % total
|
||||
end if
|
||||
|
||||
! Select smaller of the two distances
|
||||
distance = min(d_boundary, d_collision)
|
||||
|
||||
! Advance particle
|
||||
do j = 1, p % n_coord
|
||||
p % coord(j) % xyz = p % coord(j) % xyz + distance * p % coord(j) % uvw
|
||||
end do
|
||||
|
||||
! Score track-length tallies
|
||||
if (active_tracklength_tallies_size() > 0) then
|
||||
call score_tracklength_tally(p, distance)
|
||||
end if
|
||||
|
||||
! Score track-length estimate of k-eff
|
||||
if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
|
||||
global_tally_tracklength = global_tally_tracklength + p % wgt * &
|
||||
distance * material_xs % nu_fission
|
||||
end if
|
||||
|
||||
! Score flux derivative accumulators for differential tallies.
|
||||
if (active_tallies_size() > 0) call score_track_derivative(p, distance)
|
||||
|
||||
if (d_collision > d_boundary) then
|
||||
! ====================================================================
|
||||
! PARTICLE CROSSES SURFACE
|
||||
|
||||
if (next_level > 0) p % n_coord = next_level
|
||||
|
||||
! Saving previous cell data
|
||||
do j = 1, p % n_coord
|
||||
p % last_cell(j) = p % coord(j) % cell
|
||||
end do
|
||||
p % last_n_coord = p % n_coord
|
||||
|
||||
if (any(lattice_translation /= 0)) then
|
||||
! Particle crosses lattice boundary
|
||||
p % surface = ERROR_INT
|
||||
call cross_lattice(p, lattice_translation)
|
||||
p % event = EVENT_LATTICE
|
||||
else
|
||||
! Particle crosses surface
|
||||
p % surface = surface_crossed
|
||||
|
||||
call cross_surface(p)
|
||||
p % event = EVENT_SURFACE
|
||||
end if
|
||||
! Score cell to cell partial currents
|
||||
if (active_surface_tallies_size() > 0) call score_surface_tally(p)
|
||||
else
|
||||
! ====================================================================
|
||||
! PARTICLE HAS COLLISION
|
||||
|
||||
! Score collision estimate of keff
|
||||
if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then
|
||||
global_tally_collision = global_tally_collision + p % wgt * &
|
||||
material_xs % nu_fission / material_xs % total
|
||||
end if
|
||||
|
||||
! 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 (active_meshsurf_tallies_size() > 0) call score_meshsurface_tally(p)
|
||||
|
||||
! Clear surface component
|
||||
p % surface = ERROR_INT
|
||||
|
||||
if (run_CE) then
|
||||
call collision(p)
|
||||
else
|
||||
call collision_mg(p)
|
||||
end if
|
||||
|
||||
! 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 (active_collision_tallies_size() > 0) call score_collision_tally(p)
|
||||
if (active_analog_tallies_size() > 0) call score_analog_tally(p)
|
||||
|
||||
! Reset banked weight during collision
|
||||
p % n_bank = 0
|
||||
p % wgt_bank = ZERO
|
||||
p % n_delayed_bank(:) = 0
|
||||
|
||||
! Reset fission logical
|
||||
p % fission = .false.
|
||||
|
||||
! Save coordinates for tallying purposes
|
||||
p % last_xyz_current = p % coord(1) % xyz
|
||||
|
||||
! Set last material to none since cross sections will need to be
|
||||
! re-evaluated
|
||||
p % last_material = NONE
|
||||
|
||||
! Set all uvws to base level -- right now, after a collision, only the
|
||||
! base level uvws are changed
|
||||
do j = 1, p % n_coord - 1
|
||||
if (p % coord(j + 1) % rotated) then
|
||||
! If next level is rotated, apply rotation matrix
|
||||
p % coord(j + 1) % uvw = matmul(cells(p % coord(j) % cell + 1) % &
|
||||
rotation_matrix, p % coord(j) % uvw)
|
||||
else
|
||||
! Otherwise, copy this level's direction
|
||||
p % coord(j + 1) % uvw = p % coord(j) % uvw
|
||||
end if
|
||||
end do
|
||||
|
||||
! Score flux derivative accumulators for differential tallies.
|
||||
if (active_tallies_size() > 0) call score_collision_derivative(p)
|
||||
end if
|
||||
|
||||
! If particle has too many events, display warning and kill it
|
||||
n_event = n_event + 1
|
||||
if (n_event == MAX_EVENTS) then
|
||||
if (master) call warning("Particle " // trim(to_str(p%id)) &
|
||||
&// " underwent maximum number of events.")
|
||||
p % alive = .false.
|
||||
end if
|
||||
|
||||
! Check for secondary particles if this particle is dead
|
||||
if (.not. p % alive) then
|
||||
if (p % n_secondary > 0) then
|
||||
call particle_from_source(p, p % secondary_bank(p % n_secondary))
|
||||
p % n_secondary = p % n_secondary - 1
|
||||
n_event = 0
|
||||
|
||||
! Enter new particle in particle track file
|
||||
if (p % write_track) call add_particle_track()
|
||||
else
|
||||
exit EVENT_LOOP
|
||||
end if
|
||||
end if
|
||||
end do EVENT_LOOP
|
||||
|
||||
! Finish particle track output.
|
||||
if (p % write_track) then
|
||||
call write_particle_track(p)
|
||||
call finalize_particle_track(p)
|
||||
endif
|
||||
|
||||
end subroutine transport
|
||||
|
||||
!===============================================================================
|
||||
! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of
|
||||
! the geometry, is reflected, or crosses into a new lattice or cell
|
||||
!===============================================================================
|
||||
|
||||
subroutine cross_surface(p)
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
real(8) :: u ! x-component of direction
|
||||
real(8) :: v ! y-component of direction
|
||||
real(8) :: w ! z-component of direction
|
||||
real(8) :: norm ! "norm" of surface normal
|
||||
real(8) :: xyz(3) ! Saved global coordinate
|
||||
integer :: i_surface ! index in surfaces
|
||||
#ifdef DAGMC
|
||||
integer :: i_cell ! index of new cell
|
||||
#endif
|
||||
logical :: rotational ! if rotational periodic BC applied
|
||||
logical :: found ! particle found in universe?
|
||||
class(Surface), pointer :: surf
|
||||
class(Surface), pointer :: surf2 ! periodic partner surface
|
||||
|
||||
i_surface = abs(p % surface)
|
||||
surf => surfaces(i_surface)
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" Crossing surface " // trim(to_str(surf % id())))
|
||||
end if
|
||||
|
||||
if (surf % bc() == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
|
||||
! =======================================================================
|
||||
! PARTICLE LEAKS OUT OF PROBLEM
|
||||
|
||||
! Kill particle
|
||||
p % alive = .false.
|
||||
|
||||
! Score any surface current tallies -- note that the particle is moved
|
||||
! forward slightly so that if the mesh boundary is on the surface, it is
|
||||
! still processed
|
||||
|
||||
if (active_meshsurf_tallies_size() > 0) then
|
||||
! TODO: Find a better solution to score surface currents than
|
||||
! physically moving the particle forward slightly
|
||||
|
||||
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
|
||||
call score_meshsurface_tally(p)
|
||||
end if
|
||||
|
||||
! Score to global leakage tally
|
||||
global_tally_leakage = global_tally_leakage + p % wgt
|
||||
|
||||
! Display message
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" Leaked out of surface " &
|
||||
&// trim(to_str(surf % id())))
|
||||
end if
|
||||
return
|
||||
|
||||
elseif (surf % bc() == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then
|
||||
! =======================================================================
|
||||
! PARTICLE REFLECTS FROM SURFACE
|
||||
|
||||
! Do not handle reflective boundary conditions on lower universes
|
||||
if (p % n_coord /= 1) then
|
||||
call particle_mark_as_lost(p, "Cannot reflect particle " &
|
||||
// trim(to_str(p % id)) // " off surface in a lower universe.")
|
||||
return
|
||||
end if
|
||||
|
||||
! 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 (active_surface_tallies_size() > 0) call score_surface_tally(p)
|
||||
|
||||
|
||||
if (active_meshsurf_tallies_size() > 0) then
|
||||
xyz = p % coord(1) % xyz
|
||||
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
|
||||
call score_meshsurface_tally(p)
|
||||
p % coord(1) % xyz = xyz
|
||||
end if
|
||||
|
||||
! Reflect particle off surface
|
||||
call surf % reflect(p%coord(1)%xyz, p%coord(1)%uvw)
|
||||
|
||||
! Make sure new particle direction is normalized
|
||||
u = p%coord(1)%uvw(1)
|
||||
v = p%coord(1)%uvw(2)
|
||||
w = p%coord(1)%uvw(3)
|
||||
norm = sqrt(u*u + v*v + w*w)
|
||||
p%coord(1)%uvw(:) = [u, v, w] / norm
|
||||
|
||||
! Reassign particle's cell and surface
|
||||
p % coord(1) % cell = p % last_cell(p % last_n_coord)
|
||||
p % surface = -p % 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.
|
||||
|
||||
p % n_coord = 1
|
||||
call find_cell(p, found, .true.)
|
||||
if (.not. found) then
|
||||
call particle_mark_as_lost(p, "Couldn't find particle after reflecting&
|
||||
& from surface " // trim(to_str(surf % id())) // ".")
|
||||
return
|
||||
end if
|
||||
|
||||
! Set previous coordinate going slightly past surface crossing
|
||||
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
|
||||
|
||||
! Diagnostic message
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" Reflected from surface " &
|
||||
&// trim(to_str(surf%id())))
|
||||
end if
|
||||
return
|
||||
elseif (surf % bc() == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
|
||||
! =======================================================================
|
||||
! PERIODIC BOUNDARY
|
||||
|
||||
! Do not handle periodic boundary conditions on lower universes
|
||||
if (p % n_coord /= 1) then
|
||||
call particle_mark_as_lost(p, "Cannot transfer particle " &
|
||||
// trim(to_str(p % id)) // " across surface in a lower universe.&
|
||||
& Boundary conditions must be applied to universe 0.")
|
||||
return
|
||||
end if
|
||||
|
||||
! 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 (active_meshsurf_tallies_size() > 0) then
|
||||
xyz = p % coord(1) % xyz
|
||||
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
|
||||
call score_meshsurface_tally(p)
|
||||
p % coord(1) % xyz = xyz
|
||||
end if
|
||||
|
||||
! Get a pointer to the partner periodic surface. Offset the index to
|
||||
! correct for C vs. Fortran indexing.
|
||||
surf2 => surfaces(surf % i_periodic() + 1)
|
||||
|
||||
! Adjust the particle's location and direction.
|
||||
rotational = surf2 % periodic_translate(surf, p % coord(1) % xyz, &
|
||||
p % coord(1) % uvw)
|
||||
|
||||
! Reassign particle's surface
|
||||
if (rotational) then
|
||||
p % surface = surf % i_periodic() + 1
|
||||
else
|
||||
p % surface = sign(surf % i_periodic() + 1, p % surface)
|
||||
end if
|
||||
|
||||
! Figure out what cell particle is in now
|
||||
p % n_coord = 1
|
||||
call find_cell(p, found, .true.)
|
||||
if (.not. found) then
|
||||
call particle_mark_as_lost(p, "Couldn't find particle after hitting &
|
||||
&periodic boundary on surface " // trim(to_str(surf % id())) &
|
||||
// ".")
|
||||
return
|
||||
end if
|
||||
|
||||
! Set previous coordinate going slightly past surface crossing
|
||||
p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
|
||||
|
||||
! Diagnostic message
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" Hit periodic boundary on surface " &
|
||||
// trim(to_str(surf%id())))
|
||||
end if
|
||||
return
|
||||
end if
|
||||
|
||||
! ==========================================================================
|
||||
! SEARCH NEIGHBOR LISTS FOR NEXT CELL
|
||||
|
||||
#ifdef DAGMC
|
||||
if (dagmc) then
|
||||
i_cell = next_cell(cells(p % last_cell(1) + 1), surfaces(abs(p % surface)))
|
||||
! save material and temp
|
||||
p % last_material = p % material
|
||||
p % last_sqrtkT = p % sqrtKT
|
||||
! set new cell value
|
||||
p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing
|
||||
p % cell_instance = 1
|
||||
p % material = cells(i_cell) % material(1)
|
||||
p % sqrtKT = cells(i_cell) % sqrtKT(0)
|
||||
return
|
||||
end if
|
||||
#endif
|
||||
|
||||
call find_cell(p, found, .true.)
|
||||
if (found) return
|
||||
|
||||
! ==========================================================================
|
||||
! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
|
||||
|
||||
! Remove lower coordinate levels and assignment of surface
|
||||
p % surface = ERROR_INT
|
||||
p % n_coord = 1
|
||||
call find_cell(p, found)
|
||||
|
||||
if (run_mode /= MODE_PLOTTING .and. (.not. found)) then
|
||||
! 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.
|
||||
|
||||
p % n_coord = 1
|
||||
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
|
||||
call find_cell(p, found)
|
||||
|
||||
! Couldn't find next cell anywhere! This probably means there is an actual
|
||||
! undefined region in the geometry.
|
||||
|
||||
if (.not. found) then
|
||||
call particle_mark_as_lost(p, "After particle " // trim(to_str(p % id)) &
|
||||
// " crossed surface " // trim(to_str(surf % id())) &
|
||||
// " it could not be located in any cell and it did not leak.")
|
||||
return
|
||||
end if
|
||||
end if
|
||||
|
||||
end subroutine cross_surface
|
||||
|
||||
end module tracking
|
||||
Loading…
Add table
Add a link
Reference in a new issue