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

Merge pull request #1169 from paulromano/cpp-tracking

Browse source 
Translate transport, cross_surface, and particle restarts to C++
This commit is contained in:
Sterling Harper 2019-02-20 12:17:19 -05:00 • committed by GitHub
commit ef2767a127
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
61 changed files with 1001 additions and 3600 deletions
Download patch file Download diff file Expand all files Collapse all files

12
CMakeLists.txt
View file

@ -303,32 +303,23 @@ target_compile_options(faddeeva PRIVATE ${cxxflags})
#===============================================================================
add_library(libopenmc SHARED
src/algorithm.F90
src/bank_header.F90
src/api.F90
src/dagmc_header.F90
src/constants.F90
src/dict_header.F90
src/eigenvalue.F90
src/endf.F90
src/error.F90
src/geometry.F90
src/geometry_header.F90
src/hdf5_interface.F90
src/initialize.F90
src/input_xml.F90
src/material_header.F90
src/math.F90
src/mesh_header.F90
src/message_passing.F90
src/mgxs_interface.F90
src/nuclide_header.F90
src/particle_header.F90
src/particle_restart.F90
src/photon_header.F90
src/pugixml/pugixml_f.F90
src/random_lcg.F90
src/reaction_header.F90
src/relaxng
src/settings.F90
src/simulation_header.F90
@ -336,8 +327,6 @@ add_library(libopenmc SHARED
src/state_point.F90
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
@ -378,6 +367,7 @@ add_library(libopenmc SHARED
src/nuclide.cpp
src/output.cpp
src/particle.cpp
src/particle_restart.cpp
src/photon.cpp
src/physics.cpp
src/physics_common.cpp

4
include/openmc/capi.h
View file

@ -83,8 +83,7 @@ extern "C" {
int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_next_batch(int* status);
int openmc_nuclide_name(int index, char** name);
int openmc_particle_restart();
int openmc_nuclide_name(int index, const char** name);
int openmc_plot_geometry();
int openmc_reset();
int openmc_run();
@ -165,7 +164,6 @@ extern "C" {
// Global variables
extern char openmc_err_msg[256];
extern int n_nuclides;
extern int32_t n_realizations;
extern int32_t n_sab_tables;
extern int32_t n_tallies;

18
include/openmc/cell.h
View file

@ -9,14 +9,14 @@
#include "hdf5.h"
#include "pugixml.hpp"
#ifdef DAGMC
#include "DagMC.hpp"
#endif
#include "openmc/constants.h"
#include "openmc/neighbor_list.h"
#include "openmc/position.h"
#ifdef DAGMC
#include "DagMC.hpp"
#endif
#include "openmc/surface.h"
namespace openmc {
@ -195,5 +195,15 @@ public:
};
#endif
//==============================================================================
// Non-member functions
//==============================================================================
void read_cells(pugi::xml_node node);
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed);
#endif
} // namespace openmc
#endif // OPENMC_CELL_H

1
include/openmc/dagmc.h
View file

@ -26,6 +26,7 @@ extern moab::DagMC* DAG;
extern "C" void load_dagmc_geometry();
extern "C" void free_memory_dagmc();
void read_geometry_dagmc();
extern "C" pugi::xml_document* read_uwuw_materials();
bool get_uwuw_materials_xml(std::string& s);

2
include/openmc/error.h
View file

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

4
include/openmc/geometry_aux.h
View file

@ -10,6 +10,8 @@
namespace openmc {
void read_geometry_xml();
//==============================================================================
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
@ -111,7 +113,7 @@ extern "C" int maximum_levels(int32_t univ);
//! Deallocates global vectors and maps for cells, universes, and lattices.
//==============================================================================
extern "C" void free_memory_geometry_c();
extern "C" void free_memory_geometry();
} // namespace openmc
#endif // OPENMC_GEOMETRY_AUX_H

17
include/openmc/hdf5_interface.h
View file

@ -223,8 +223,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
// Templates/overloads for read_dataset and related methods
//==============================================================================
template<typename T>
void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
// Template for scalars. We need to be careful that the compiler does not use
// this version of read_dataset for vectors, arrays, or other non-scalar types.
// enable_if_t allows us to conditionally remove the function from overload
// resolution when the type T doesn't meet a certain criterion.
template<typename T> inline
std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
{
read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
}
@ -242,6 +247,14 @@ read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false)
str = std::string{buffer, n};
}
// array version
template<typename T, std::size_t N> inline void
read_dataset(hid_t dset, const char* name, std::array<T, N>& buffer, bool indep=false)
{
read_dataset(dset, name, H5TypeMap<T>::type_id, buffer.data(), indep);
}
// vector version
template <typename T>
void read_dataset(hid_t dset, std::vector<T>& vec, bool indep=false)
{

6
include/openmc/lattice.h
View file

@ -273,5 +273,11 @@ private:
std::array<double, 2> pitch_; //!< Lattice tile width and height
};
//==============================================================================
// Non-member functions
//==============================================================================
void read_lattices(pugi::xml_node node);
} // namespace openmc
#endif // OPENMC_LATTICE_H

6
include/openmc/particle.h
View file

@ -148,6 +148,12 @@ extern "C" {
//! \param src Source site data
void from_source(const Bank* src);
//! Transport a particle from birth to death
void transport();
//! Cross a surface and handle boundary conditions
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);

10
include/openmc/particle_restart.h Normal file
View file

@ -0,0 +1,10 @@
#ifndef OPENMC_PARTICLE_RESTART_H
#define OPENMC_PARTICLE_RESTART_H
namespace openmc {
void run_particle_restart();
} // namespace openmc
#endif // OPENMC_PARTICLE_RESTART_H

2
include/openmc/physics.h
View file

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

2
include/openmc/physics_mg.h
View file

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

4
include/openmc/simulation.h
View file

@ -36,6 +36,7 @@ extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" int total_gen; //!< total number of generations simulated
extern "C" double total_weight; //!< Total source weight in a batch
extern "C" int64_t work; //!< number of particles per process
extern std::vector<double> k_generation;
@ -83,6 +84,9 @@ void finalize_generation();
//! Determine overall generation number
extern "C" int overall_generation();
extern "C" void simulation_init_f();
extern "C" void simulation_finalize_f();
#ifdef OPENMC_MPI
void broadcast_results();
#endif

6
include/openmc/surface.h
View file

@ -395,6 +395,12 @@ public:
void to_hdf5_inner(hid_t group_id) const;
};
//==============================================================================
// Non-member functions
//==============================================================================
void read_surfaces(pugi::xml_node node);
//==============================================================================
// Fortran compatibility functions
//==============================================================================

19
include/openmc/tallies/derivative.h
View file

@ -35,6 +35,25 @@ void
apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
double atom_density, int score_bin, double& score);
//! Adjust diff tally flux derivatives for a particle scattering event.
//
//! Note that this subroutine will be called after absorption events in
//! addition to scattering events, but any flux derivatives scored after an
//! absorption will never be tallied. The paricle will be killed before any
//! further tallies are scored.
//
//! \param p The particle being tracked
void score_collision_derivative(const Particle* p);
//! Adjust diff tally flux derivatives for a particle tracking event.
//
//! \param p The particle being tracked
//! \param distance The distance in [cm] traveled by the particle
void score_track_derivative(const Particle* p, double distance);
//! Set the flux derivatives on differential tallies to zero.
void zero_flux_derivs();
} // namespace openmc
//==============================================================================

2
include/openmc/tallies/tally.h
View file

@ -101,8 +101,6 @@ private:
// Global variable declarations
//==============================================================================
extern "C" double total_weight;
namespace model {
extern std::vector<std::unique_ptr<Tally>> tallies;

41
include/openmc/tallies/tally_scoring.h
View file

@ -50,6 +50,47 @@ private:
// Non-member functions
//==============================================================================
//! Score tallies using a 1 / Sigma_t estimate of the flux.
//
//! This is triggered after every collision. It is invalid for tallies that
//! require post-collison information because it can score reactions that didn't
//! actually occur, and we don't a priori know what the outcome will be for
//! reactions that we didn't sample. It is assumed the material is not void
//! since collisions do not occur in voids.
//
//! \param p The particle being tracked
void score_collision_tally(const Particle* p);
//! Score tallies based on a simple count of events (for continuous energy).
//
//! Analog tallies are triggered at every collision, not every event.
//
//! \param p The particle being tracked
void score_analog_tally_ce(const Particle* p);
//! Score tallies based on a simple count of events (for multigroup).
//
//! Analog tallies are triggered at every collision, not every event.
//
//! \param p The particle being tracked
void score_analog_tally_mg(const Particle* p);
//! Score tallies using a tracklength estimate of the flux.
//
//! This is triggered at every event (surface crossing, lattice crossing, or
//! collision) and thus cannot be done for tallies that require post-collision
//! information.
//
//! \param p The particle being tracked
//! \param distance The distance in [cm] traveled by the particle
void score_tracklength_tally(const Particle* p, double distance);
//! Score surface or mesh-surface tallies for particle currents.
//
//! \param p The particle being tracked
//! \param tallies A vector of tallies to score to
void score_surface_tally(const Particle* p, const std::vector<int>& tallies);
} // namespace openmc
#endif // OPENMC_TALLIES_TALLY_SCORING_H

5
openmc/capi/cell.py
View file

@ -38,6 +38,7 @@ _dll.openmc_cell_set_temperature.errcheck = _error_handler
_dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_cell_index.restype = c_int
_dll.openmc_get_cell_index.errcheck = _error_handler
_dll.cells_size.restype = c_int
class Cell(_FortranObjectWithID):
@ -156,10 +157,10 @@ class _CellMapping(Mapping):
def __iter__(self):
for i in range(len(self)):
yield Cell(index=i + 1).id
yield Cell(index=i).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_cells').value
return _dll.cells_size()
def __repr__(self):
return repr(dict(self))

2
openmc/capi/material.py
View file

@ -163,7 +163,7 @@ class Material(_FortranObjectWithID):
_dll.openmc_material_get_densities(self._index, nuclides, densities, n)
# Convert to appropriate types and return
nuclide_list = [Nuclide(nuclides[i] + 1).name for i in range(n.value)]
nuclide_list = [Nuclide(nuclides[i]).name for i in range(n.value)]
density_array = as_array(densities, (n.value,))
return nuclide_list, density_array

12
openmc/capi/nuclide.py
View file

@ -23,6 +23,7 @@ _dll.openmc_load_nuclide.errcheck = _error_handler
_dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)]
_dll.openmc_nuclide_name.restype = c_int
_dll.openmc_nuclide_name.errcheck = _error_handler
_dll.nuclides_size.restype = c_int
def load_nuclide(name):
@ -70,12 +71,7 @@ class Nuclide(_FortranObject):
def name(self):
name = c_char_p()
_dll.openmc_nuclide_name(self._index, name)
# Find blank in name
i = 0
while name.value[i:i+1] != b' ':
i += 1
return name.value[:i].decode()
return name.value.decode()
class _NuclideMapping(Mapping):
@ -91,10 +87,10 @@ class _NuclideMapping(Mapping):
def __iter__(self):
for i in range(len(self)):
yield Nuclide(i + 1).name
yield Nuclide(i).name
def __len__(self):
return c_int.in_dll(_dll, 'n_nuclides').value
return _dll.nuclides_size()
def __repr__(self):
return repr(dict(self))

2
openmc/capi/tally.py
View file

@ -278,7 +278,7 @@ class Tally(_FortranObjectWithID):
nucs = POINTER(c_int)()
n = c_int()
_dll.openmc_tally_get_nuclides(self._index, nucs, n)
return [Nuclide(nucs[i]+1).name if nucs[i] >= 0 else 'total'
return [Nuclide(nucs[i]).name if nucs[i] >= 0 else 'total'
for i in range(n.value)]
@nuclides.setter

141
src/algorithm.F90
View file

@ -1,141 +0,0 @@
module algorithm
use constants
use stl_vector, only: VectorInt, VectorReal
implicit none
interface sort
module procedure sort_int, sort_real, sort_vector_int, sort_vector_real
end interface sort
interface find
module procedure find_int, find_real, find_vector_int, find_vector_real
end interface find
contains
!===============================================================================
! SORT sorts an array in place using an insertion sort.
!===============================================================================
pure subroutine sort_int(array)
integer, intent(inout) :: array(:)
integer :: k, m
integer :: temp
if (size(array) > 1) then
SORT: do k = 2, size(array)
! Save value to move
m = k
temp = array(k)
MOVE_OVER: do while (m > 1)
! Check if insertion value is greater than (m-1)th value
if (temp >= array(m - 1)) exit
! Move values over until hitting one that's not larger
array(m) = array(m - 1)
m = m - 1
end do MOVE_OVER
! Put the original value into its new position
array(m) = temp
end do SORT
end if
end subroutine sort_int
pure subroutine sort_real(array)
real(8), intent(inout) :: array(:)
integer :: k, m
real(8) :: temp
if (size(array) > 1) then
SORT: do k = 2, size(array)
! Save value to move
m = k
temp = array(k)
MOVE_OVER: do while (m > 1)
! Check if insertion value is greater than (m-1)th value
if (temp >= array(m - 1)) exit
! Move values over until hitting one that's not larger
array(m) = array(m - 1)
m = m - 1
end do MOVE_OVER
! Put the original value into its new position
array(m) = temp
end do SORT
end if
end subroutine sort_real
pure subroutine sort_vector_int(vec)
type(VectorInt), intent(inout) :: vec
call sort_int(vec % data(1:vec%size()))
end subroutine sort_vector_int
pure subroutine sort_vector_real(vec)
type(VectorReal), intent(inout) :: vec
call sort_real(vec % data(1:vec%size()))
end subroutine sort_vector_real
!===============================================================================
! FIND determines the index of the first occurrence of a value in an array. If
! the value does not appear in the array, -1 is returned.
!===============================================================================
pure function find_int(array, val) result(index)
integer, intent(in) :: array(:)
integer, intent(in) :: val
integer :: index
integer :: i
index = -1
do i = 1, size(array)
if (array(i) == val) then
index = i
exit
end if
end do
end function find_int
pure function find_real(array, val) result(index)
real(8), intent(in) :: array(:)
real(8), intent(in) :: val
integer :: index
integer :: i
index = -1
do i = 1, size(array)
if (array(i) == val) then
index = i
exit
end if
end do
end function find_real
pure function find_vector_int(vec, val) result(index)
type(VectorInt), intent(in) :: vec
integer, intent(in) :: val
integer :: index
index = find_int(vec % data(1:vec % size()), val)
end function find_vector_int
pure function find_vector_real(vec, val) result(index)
type(VectorReal), intent(in) :: vec
real(8), intent(in) :: val
integer :: index
index = find_real(vec % data(1:vec % size()), val)
end function find_vector_real
end module algorithm

11
src/api.F90
View file

@ -39,7 +39,7 @@ contains
err = E_UNASSIGNED
if (found) then
index = p % coord(p % n_coord) % cell + 1
index = p % coord(p % n_coord) % cell
instance = p % cell_instance
err = 0
else
@ -59,12 +59,11 @@ contains
subroutine free_memory() bind(C)
use bank_header
use geometry_header
use material_header
use nuclide_header
use photon_header
use settings
use simulation_header
use surface_header
use tally_filter_header
use tally_header
@ -90,6 +89,12 @@ contains
subroutine free_memory_volume() bind(C)
end subroutine
subroutine free_memory_surfaces() bind(C)
end subroutine
subroutine free_memory_geometry() bind(C)
end subroutine
subroutine sab_clear() bind(C)
end subroutine
end interface

160
src/cell.cpp
View file

@ -23,8 +23,6 @@ namespace openmc {
namespace model {
int32_t n_cells {0};
std::vector<Cell*> cells;
std::unordered_map<int32_t, int32_t> cell_map;
@ -643,18 +641,18 @@ void DAGCell::to_hdf5(hid_t group_id) const { return; }
// Non-method functions
//==============================================================================
extern "C" void
read_cells(pugi::xml_node* node)
void read_cells(pugi::xml_node node)
{
// Count the number of cells.
for (pugi::xml_node cell_node: node->children("cell")) {model::n_cells++;}
if (model::n_cells == 0) {
int n_cells = 0;
for (pugi::xml_node cell_node: node.children("cell")) {n_cells++;}
if (n_cells == 0) {
fatal_error("No cells found in geometry.xml!");
}
// Loop over XML cell elements and populate the array.
model::cells.reserve(model::n_cells);
for (pugi::xml_node cell_node : node->children("cell")) {
model::cells.reserve(n_cells);
for (pugi::xml_node cell_node : node.children("cell")) {
model::cells.push_back(new CSGCell(cell_node));
}
@ -699,9 +697,8 @@ read_cells(pugi::xml_node* node)
extern "C" int
openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
{
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
*type = c.type_;
if (c.type_ == FILL_MATERIAL) {
*indices = c.material_.data();
@ -721,9 +718,8 @@ extern "C" int
openmc_cell_set_fill(int32_t index, int type, int32_t n,
const int32_t* indices)
{
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
if (type == FILL_MATERIAL) {
c.type_ = FILL_MATERIAL;
c.material_.clear();
@ -756,9 +752,8 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
extern "C" int
openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
{
if (index >= 1 && index <= model::cells.size()) {
//TODO: off-by-one
Cell& c {*model::cells[index - 1]};
if (index >= 0 && index < model::cells.size()) {
Cell& c {*model::cells[index]};
if (instance) {
if (*instance >= 0 && *instance < c.sqrtkT_.size()) {
@ -781,71 +776,78 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
return 0;
}
//! Return the index in the cells array of a cell with a given ID
extern "C" int
openmc_get_cell_index(int32_t id, int32_t* index)
{
auto it = model::cell_map.find(id);
if (it != model::cell_map.end()) {
*index = it->second;
return 0;
} else {
set_errmsg("No cell exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
}
}
//! Return the ID of a cell
extern "C" int
openmc_cell_get_id(int32_t index, int32_t* id)
{
if (index >= 0 && index < model::cells.size()) {
*id = model::cells[index]->id_;
return 0;
} else {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//! Set the ID of a cell
extern "C" int
openmc_cell_set_id(int32_t index, int32_t id)
{
if (index >= 0 && index < model::cells.size()) {
model::cells[index]->id_ = id;
model::cell_map[id] = index;
return 0;
} else {
set_errmsg("Index in cells array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//! Extend the cells array by n elements
extern "C" int
openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end)
{
if (index_start) *index_start = model::cells.size();
if (index_end) *index_end = model::cells.size() + n - 1;
for (int32_t i = 0; i < n; i++) {
model::cells.push_back(new CSGCell());
}
return 0;
}
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed)
{
moab::EntityHandle surf =
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_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol);
}
#endif
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Cell* cell_pointer(int32_t cell_ind) {return model::cells[cell_ind];}
int32_t cell_id(Cell* c) {return c->id_;}
void
cell_set_id(Cell* c, int32_t id)
{
c->id_ = id;
// Find the index of this cell and update the cell map.
for (int i = 0; i < model::cells.size(); i++) {
if (model::cells[i] == c) {
model::cell_map[id] = i;
break;
}
}
}
#ifdef DAGMC
int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed )
{
moab::EntityHandle surf =
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_);
moab::EntityHandle new_vol;
cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol);
return cur_cell->dagmc_ptr_->index_by_handle(new_vol);
}
#endif
int32_t cell_universe(Cell* c) {return c->universe_;}
int32_t cell_fill(Cell* c) {return c->fill_;}
int cell_material_size(Cell* c) {return c->material_.size();}
//TODO: off-by-one
int32_t cell_material(Cell* c, int i)
{
int32_t mat = c->material_[i-1];
if (mat == MATERIAL_VOID) return MATERIAL_VOID;
return mat + 1;
}
int cell_sqrtkT_size(Cell* c) {return c->sqrtkT_.size();}
double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];}
void extend_cells_c(int32_t n)
{
model::cells.reserve(model::cells.size() + n);
for (int32_t i = 0; i < n; i++) {
model::cells.push_back(new CSGCell());
}
model::n_cells = model::cells.size();
}
}
extern "C" int cells_size() { return model::cells.size(); }
} // namespace openmc

121
src/constants.F90
View file

@ -7,16 +7,8 @@ module constants
! ============================================================================
! VERSIONING NUMBERS
! OpenMC major, minor, and release numbers
integer, parameter :: VERSION_MAJOR = 0
integer, parameter :: VERSION_MINOR = 10
integer, parameter :: VERSION_RELEASE = 0
integer, parameter :: &
VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE]
! Version numbers for binary files
integer, parameter :: VERSION_TRACK(2) = [2, 0]
integer, parameter :: VERSION_VOLUME(2) = [1, 0]
integer, parameter :: VERSION_TRACK(2) = [2, 0]
! ============================================================================
! ADJUSTABLE PARAMETERS
@ -24,15 +16,6 @@ module constants
! NOTE: This is the only section of the constants module that should ever be
! adjusted. Modifying constants in other sections may cause the code to fail.
! Significance level for confidence intervals
real(8), parameter :: CONFIDENCE_LEVEL = 0.95_8
! Used for surface current tallies
real(8), parameter :: TINY_BIT = 1e-8_8
! Maximum number of collisions/crossings
integer, parameter :: MAX_EVENTS = 1000000
! Maximum number of secondary particles created
integer, parameter :: MAX_SECONDARY = 1000
@ -50,24 +33,12 @@ module constants
real(8), parameter :: &
PI = 3.1415926535898_8, & ! pi
MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu
MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2
C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s
K_BOLTZMANN = 8.6173303e-5_8, & ! Boltzmann constant in eV/K
INFINITY = huge(0.0_8), & ! positive infinity
ZERO = 0.0_8, &
HALF = 0.5_8, &
ONE = 1.0_8, &
TWO = 2.0_8
! Electron subshell labels
character(3), parameter :: SUBSHELLS(39) = [ &
'K ', 'L1 ', 'L2 ', 'L3 ', 'M1 ', 'M2 ', 'M3 ', 'M4 ', 'M5 ', &
'N1 ', 'N2 ', 'N3 ', 'N4 ', 'N5 ', 'N6 ', 'N7 ', 'O1 ', 'O2 ', &
'O3 ', 'O4 ', 'O5 ', 'O6 ', 'O7 ', 'O8 ', 'O9 ', 'P1 ', 'P2 ', &
'P3 ', 'P4 ', 'P5 ', 'P6 ', 'P7 ', 'P8 ', 'P9 ', 'P10', 'P11', &
'Q1 ', 'Q2 ', 'Q3 ']
! ============================================================================
! GEOMETRY-RELATED CONSTANTS
@ -99,55 +70,11 @@ module constants
ELECTRON = 2, &
POSITRON = 3
! Reaction types
integer, parameter :: &
TOTAL_XS = 1, ELASTIC = 2, N_NONELASTIC = 3, N_LEVEL = 4, MISC = 5, &
N_2ND = 11, N_2N = 16, N_3N = 17, N_FISSION = 18, N_F = 19, &
N_NF = 20, N_2NF = 21, N_NA = 22, N_N3A = 23, N_2NA = 24, &
N_3NA = 25, N_NP = 28, N_N2A = 29, N_2N2A = 30, N_ND = 32, &
N_NT = 33, N_N3HE = 34, N_ND2A = 35, N_NT2A = 36, N_4N = 37, &
N_3NF = 38, N_2NP = 41, N_3NP = 42, N_N2P = 44, N_NPA = 45, &
N_N1 = 51, N_N40 = 90, N_NC = 91, N_DISAPPEAR = 101, N_GAMMA = 102, &
N_P = 103, N_D = 104, N_T = 105, N_3HE = 106, N_A = 107, &
N_2A = 108, N_3A = 109, N_2P = 111, N_PA = 112, N_T2A = 113, &
N_D2A = 114, N_PD = 115, N_PT = 116, N_DA = 117, N_5N = 152, &
N_6N = 153, N_2NT = 154, N_TA = 155, N_4NP = 156, N_3ND = 157, &
N_NDA = 158, N_2NPA = 159, N_7N = 160, N_8N = 161, N_5NP = 162, &
N_6NP = 163, N_7NP = 164, N_4NA = 165, N_5NA = 166, N_6NA = 167, &
N_7NA = 168, N_4ND = 169, N_5ND = 170, N_6ND = 171, N_3NT = 172, &
N_4NT = 173, N_5NT = 174, N_6NT = 175, N_2N3HE = 176, N_3N3HE = 177, &
N_4N3HE = 178, N_3N2P = 179, N_3N3A = 180, N_3NPA = 181, N_DT = 182, &
N_NPD = 183, N_NPT = 184, N_NDT = 185, N_NP3HE = 186, N_ND3HE = 187, &
N_NT3HE = 188, N_NTA = 189, N_2N2P = 190, N_P3HE = 191, N_D3HE = 192, &
N_3HEA = 193, N_4N2P = 194, N_4N2A = 195, N_4NPA = 196, N_3P = 197, &
N_N3P = 198, N_3N2PA = 199, N_5N2P = 200, N_P0 = 600, N_PC = 649, &
N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, &
N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891, &
COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, &
PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517
! Depletion reactions
integer, parameter :: DEPLETION_RX(6) = [N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N]
! MGXS Table Types
integer, parameter :: &
MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data
MGXS_ANGLE = 2 ! Data by Angular Bins
! Secondary particle emission type
integer, parameter :: &
EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle
EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle
EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed)
! Maximum number of partial fission reactions
integer, parameter :: PARTIAL_FISSION_MAX = 4
! Temperature treatment method
integer, parameter :: &
TEMPERATURE_NEAREST = 1, &
TEMPERATURE_INTERPOLATION = 2
! ============================================================================
! TALLY-RELATED CONSTANTS
@ -206,21 +133,6 @@ module constants
FILTER_MATERIAL = 2, &
FILTER_CELL = 3
! Tally surface current directions
integer, parameter :: &
OUT_LEFT = 1, & ! x min
IN_LEFT = 2, & ! x min
OUT_RIGHT = 3, & ! x max
IN_RIGHT = 4, & ! x max
OUT_BACK = 5, & ! y min
IN_BACK = 6, & ! y min
OUT_FRONT = 7, & ! y max
IN_FRONT = 8, & ! y max
OUT_BOTTOM = 9, & ! z min
IN_BOTTOM = 10, & ! z min
OUT_TOP = 11, & ! z max
IN_TOP = 12 ! z max
! Global tally parameters
integer, parameter :: N_GLOBAL_TALLIES = 4
integer, parameter :: &
@ -235,37 +147,6 @@ module constants
DIFF_NUCLIDE_DENSITY = 2, &
DIFF_TEMPERATURE = 3
! Mgxs::get_xs enumerated types
integer(C_INT), parameter :: &
MG_GET_XS_TOTAL = 0, &
MG_GET_XS_ABSORPTION = 1, &
MG_GET_XS_INVERSE_VELOCITY = 2, &
MG_GET_XS_DECAY_RATE = 3, &
MG_GET_XS_SCATTER = 4, &
MG_GET_XS_SCATTER_MULT = 5, &
MG_GET_XS_SCATTER_FMU_MULT = 6, &
MG_GET_XS_SCATTER_FMU = 7, &
MG_GET_XS_FISSION = 8, &
MG_GET_XS_KAPPA_FISSION = 9, &
MG_GET_XS_PROMPT_NU_FISSION = 10, &
MG_GET_XS_DELAYED_NU_FISSION = 11, &
MG_GET_XS_NU_FISSION = 12, &
MG_GET_XS_CHI_PROMPT = 13, &
MG_GET_XS_CHI_DELAYED = 14
! ============================================================================
! RANDOM NUMBER STREAM CONSTANTS
integer(C_INT), bind(C, name='N_STREAMS') :: N_STREAMS
integer(C_INT), bind(C, name='STREAM_TRACKING') :: STREAM_TRACKING
integer(C_INT), bind(C, name='STREAM_TALLIES') :: STREAM_TALLIES
integer(C_INT), bind(C, name='STREAM_SOURCE') :: STREAM_SOURCE
integer(C_INT), bind(C, name='STREAM_URR_PTABLE') :: STREAM_URR_PTABLE
integer(C_INT), bind(C, name='STREAM_VOLUME') :: STREAM_VOLUME
integer(C_INT), bind(C, name='STREAM_PHOTON') :: STREAM_PHOTON
integer(C_INT64_T), parameter :: DEFAULT_SEED = 1_8
! ============================================================================
! MISCELLANEOUS CONSTANTS

7
src/cross_sections.cpp
View file

@ -179,9 +179,6 @@ void read_cross_sections_xml()
}
}
extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr,
const double* temps, int n, int n_nuclide);
void
read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
const std::vector<std::vector<double>>& thermal_temps)
@ -226,10 +223,6 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
data::nuclides.push_back(std::make_unique<Nuclide>(
group, nuc_temps[i_nuc], i_nuclide));
// Read from Fortran too
nuclide_from_hdf5(group, data::nuclides.back().get(),
&nuc_temps[i_nuc].front(), nuc_temps[i_nuc].size(), i_nuclide + 1);
close_group(group);
file_close(file_id);

19
src/dagmc.cpp
View file

@ -4,6 +4,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry_aux.h"
#include "openmc/string_utils.h"
#include "openmc/settings.h"
#include "openmc/geometry.h"
@ -133,9 +134,9 @@ void load_dagmc_geometry()
/// Cells (Volumes) \\\
// initialize cell objects
model::n_cells = model::DAG->num_entities(3);
int n_cells = model::DAG->num_entities(3);
moab::EntityHandle graveyard = 0;
for (int i = 0; i < model::n_cells; i++) {
for (int i = 0; i < n_cells; i++) {
moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1);
// set cell ids using global IDs
@ -312,6 +313,20 @@ void load_dagmc_geometry()
return;
}
void read_geometry_dagmc()
{
// Check if dagmc.h5m exists
std::string filename = settings::path_input + "dagmc.h5m";
if (!file_exists(filename)) {
fatal_error("Geometry DAGMC file '" + filename + "' does not exist!");
}
write_message("Reading DAGMC geometry...", 5);
load_dagmc_geometry();
model::root_universe = find_root_universe();
}
void free_memory_dagmc()
{
delete model::DAG;

15
src/eigenvalue.F90
View file

@ -1,15 +0,0 @@
module eigenvalue
use, intrinsic :: ISO_C_BINDING
implicit none
interface
function openmc_get_keff(k_combined) result(err) bind(C)
import C_INT, C_DOUBLE
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
end function
end interface
end module eigenvalue

268
src/endf.F90
View file

@ -1,268 +0,0 @@
module endf
use constants
use string, only: to_str
implicit none
contains
!===============================================================================
! REACTION_NAME gives the name of the reaction for a given MT value
!===============================================================================
pure function reaction_name(MT) result(string)
integer, intent(in) :: MT
character(MAX_WORD_LEN) :: string
select case (MT)
! Special reactions for tallies
case (SCORE_FLUX)
string = "flux"
case (SCORE_TOTAL)
string = "total"
case (SCORE_SCATTER)
string = "scatter"
case (SCORE_NU_SCATTER)
string = "nu-scatter"
case (SCORE_ABSORPTION)
string = "absorption"
case (SCORE_FISSION)
string = "fission"
case (SCORE_NU_FISSION)
string = "nu-fission"
case (SCORE_DECAY_RATE)
string = "decay-rate"
case (SCORE_DELAYED_NU_FISSION)
string = "delayed-nu-fission"
case (SCORE_PROMPT_NU_FISSION)
string = "prompt-nu-fission"
case (SCORE_KAPPA_FISSION)
string = "kappa-fission"
case (SCORE_CURRENT)
string = "current"
case (SCORE_EVENTS)
string = "events"
case (SCORE_INVERSE_VELOCITY)
string = "inverse-velocity"
case (SCORE_FISS_Q_PROMPT)
string = "fission-q-prompt"
case (SCORE_FISS_Q_RECOV)
string = "fission-q-recoverable"
! Normal ENDF-based reactions
case (TOTAL_XS)
string = '(n,total)'
case (ELASTIC)
string = '(n,elastic)'
case (N_LEVEL)
string = '(n,level)'
case (N_2ND)
string = '(n,2nd)'
case (N_2N)
string = '(n,2n)'
case (N_3N)
string = '(n,3n)'
case (N_FISSION)
string = '(n,fission)'
case (N_F)
string = '(n,f)'
case (N_NF)
string = '(n,nf)'
case (N_2NF)
string = '(n,2nf)'
case (N_NA)
string = '(n,na)'
case (N_N3A)
string = '(n,n3a)'
case (N_2NA)
string = '(n,2na)'
case (N_3NA)
string = '(n,3na)'
case (N_NP)
string = '(n,np)'
case (N_N2A)
string = '(n,n2a)'
case (N_2N2A)
string = '(n,2n2a)'
case (N_ND)
string = '(n,nd)'
case (N_NT)
string = '(n,nt)'
case (N_N3HE)
string = '(n,nHe-3)'
case (N_ND2A)
string = '(n,nd2a)'
case (N_NT2A)
string = '(n,nt2a)'
case (N_4N)
string = '(n,4n)'
case (N_3NF)
string = '(n,3nf)'
case (N_2NP)
string = '(n,2np)'
case (N_3NP)
string = '(n,3np)'
case (N_N2P)
string = '(n,n2p)'
case (N_NPA)
string = '(n,npa)'
case (N_N1 : N_N40)
string = '(n,n' // trim(to_str(MT-50)) // ')'
case (N_NC)
string = '(n,nc)'
case (N_DISAPPEAR)
string = '(n,disappear)'
case (N_GAMMA)
string = '(n,gamma)'
case (N_P)
string = '(n,p)'
case (N_D)
string = '(n,d)'
case (N_T)
string = '(n,t)'
case (N_3HE)
string = '(n,3He)'
case (N_A)
string = '(n,a)'
case (N_2A)
string = '(n,2a)'
case (N_3A)
string = '(n,3a)'
case (N_2P)
string = '(n,2p)'
case (N_PA)
string = '(n,pa)'
case (N_T2A)
string = '(n,t2a)'
case (N_D2A)
string = '(n,d2a)'
case (N_PD)
string = '(n,pd)'
case (N_PT)
string = '(n,pt)'
case (N_DA)
string = '(n,da)'
case (201)
string = '(n,Xn)'
case (202)
string = '(n,Xgamma)'
case (203)
string = '(n,Xp)'
case (204)
string = '(n,Xd)'
case (205)
string = '(n,Xt)'
case (206)
string = '(n,X3He)'
case (207)
string = '(n,Xa)'
case (444)
string = '(damage)'
case (COHERENT)
string = 'coherent scatter'
case (INCOHERENT)
string = 'incoherent scatter'
case (PAIR_PROD_ELEC)
string = 'pair production, electron'
case (PAIR_PROD)
string = 'pair production'
case (PAIR_PROD_NUC)
string = 'pair production, nuclear'
case (PHOTOELECTRIC)
string = 'photoelectric'
case (534 : 572)
string = 'photoelectric, ' // trim(SUBSHELLS(MT - 533)) // ' subshell'
case (600 : 648)
string = '(n,p' // trim(to_str(MT-600)) // ')'
case (649)
string = '(n,pc)'
case (650 : 698)
string = '(n,d' // trim(to_str(MT-650)) // ')'
case (699)
string = '(n,dc)'
case (700 : 748)
string = '(n,t' // trim(to_str(MT-700)) // ')'
case (749)
string = '(n,tc)'
case (750 : 798)
string = '(n,3He' // trim(to_str(MT-750)) // ')'
case (799)
string = '(n,3Hec)'
case (800 : 848)
string = '(n,a' // trim(to_str(MT-800)) // ')'
case (849)
string = '(n,ac)'
case default
string = 'MT=' // trim(to_str(MT))
end select
end function reaction_name
!===============================================================================
! IS_FISSION determines if a given MT number is that of a fission event. This
! accounts for aggregate fission (MT=18) as well as partial fission reactions.
!===============================================================================
function is_fission(MT) result(fission_event)
integer, intent(in) :: MT
logical :: fission_event
if (MT == N_FISSION .or. MT == N_F .or. MT == N_NF .or. MT == N_2NF &
.or. MT == N_3NF) then
fission_event = .true.
else
fission_event = .false.
end if
end function is_fission
!===============================================================================
! IS_DISAPPEARANCE determines if a given MT number is that of a disappearance
! reaction, i.e. a reaction with no neutron in the exit channel
!===============================================================================
function is_disappearance(MT) result(dis)
integer, intent(in) :: MT
logical :: dis
if (MT >= N_DISAPPEAR .and. MT <= N_DA) then
dis = .true.
elseif (MT >= N_P0 .and. MT <= N_AC) then
dis = .true.
elseif (any(MT == [N_TA, N_DT, N_P3HE, N_D3HE, N_3HEA, N_3P])) then
dis = .true.
else
dis = .false.
end if
end function is_disappearance
!===============================================================================
! IS_INELASTIC_SCATTER determines if a given MT number is that of an inelastic
! scattering event
!===============================================================================
function is_inelastic_scatter(MT) result(retval)
integer, intent(in) :: MT
logical :: retval
if (MT < 100) then
if (is_fission(MT)) then
retval = .false.
else
retval = (MT >= MISC .and. MT /= 27)
end if
elseif (MT <= 200) then
retval = (.not. is_disappearance(MT))
else
retval = .false.
end if
end function
end module endf

4
src/error.cpp
View file

@ -33,7 +33,7 @@ void output(const std::string& message, std::ostream& out, int indent)
while (i_start < length) {
if (length - i_start < line_len) {
// Remainder of message will fit on line
std::cerr << message.substr(i_start) << '\n';
out << message.substr(i_start) << '\n';
break;
} else {
@ -42,7 +42,7 @@ void output(const std::string& message, std::ostream& out, int indent)
auto pos = s.find_last_of(' ');
// Write up to last space, or whole line if no space is present
std::cerr << s.substr(0, pos) << '\n' << std::setw(indent) << " ";
out << s.substr(0, pos) << '\n' << std::setw(indent) << " ";
// Advance starting position
i_start += (pos == std::string::npos) ? line_len : pos + 1;

27
src/geometry.F90
View file

@ -1,12 +1,9 @@
module geometry
use constants
use geometry_header
use particle_header
use simulation_header
use settings
use surface_header
use string, only: to_str
use, intrinsic :: ISO_C_BINDING
@ -50,34 +47,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

44
src/geometry_aux.cpp
View file

@ -4,10 +4,14 @@
#include <sstream>
#include <unordered_set>
#include "pugixml.hpp"
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/dagmc.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/geometry.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
@ -19,6 +23,43 @@
namespace openmc {
void read_geometry_xml()
{
#ifdef DAGMC
if (settings::dagmc) {
read_geometry_dagmc();
return;
}
#endif
// Display output message
write_message("Reading geometry XML file...", 5);
// Check if geometry.xml exists
std::string filename = settings::path_input + "geometry.xml";
if (!file_exists(filename)) {
fatal_error("Geometry XML file '" + filename + "' does not exist!");
}
// Parse settings.xml file
pugi::xml_document doc;
auto result = doc.load_file(filename.c_str());
if (!result) {
fatal_error("Error processing geometry.xml file.");
}
// Get root element
pugi::xml_node root = doc.document_element();
// Read surfaces, cells, lattice
read_surfaces(root);
read_cells(root);
read_lattices(root);
// Allocate universes, universe cell arrays, and assign base universe
model::root_universe = find_root_universe();
}
//==============================================================================
void
@ -473,12 +514,11 @@ maximum_levels(int32_t univ)
//==============================================================================
void
free_memory_geometry_c()
free_memory_geometry()
{
for (Cell* c : model::cells) {delete c;}
model::cells.clear();
model::cell_map.clear();
model::n_cells = 0;
for (Universe* u : model::universes) {delete u;}
model::universes.clear();

283
src/geometry_header.F90
View file

@ -1,283 +0,0 @@
module geometry_header
use, intrinsic :: ISO_C_BINDING
use dict_header, only: DictIntInt
use nuclide_header
implicit none
interface
function cell_pointer(cell_ind) bind(C) result(ptr)
import C_PTR, C_INT32_T
integer(C_INT32_T), intent(in), value :: cell_ind
type(C_PTR) :: ptr
end function cell_pointer
function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: id
end function cell_id_c
subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id')
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T), intent(in), value :: id
end subroutine cell_set_id_c
function cell_universe_c(cell_ptr) bind(C, name='cell_universe') &
result(universe)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: universe
end function cell_universe_c
function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill)
import C_PTR, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT32_T) :: fill
end function cell_fill_c
function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') &
result(n)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT) :: n
end function cell_material_size_c
function cell_material_c(cell_ptr, i) bind(C, name='cell_material') &
result(mat)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: i
integer(C_INT32_T) :: mat
end function cell_material_c
function cell_sqrtkT_size_c(cell_ptr) bind(C, name='cell_sqrtkT_size') &
result(n)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT) :: n
end function cell_sqrtkT_size_c
function cell_sqrtkT_c(cell_ptr, i) bind(C, name='cell_sqrtkT') &
result(sqrtkT)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: i
real(C_DOUBLE) :: sqrtkT
end function cell_sqrtkT_c
subroutine extend_cells_c(n) bind(C)
import C_INT32_t
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_cells_c
end interface
!===============================================================================
! CELL defines a closed volume by its bounding surfaces
!===============================================================================
type Cell
type(C_PTR) :: ptr
! Rotation matrix and translation vector
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
contains
procedure :: id => cell_id
procedure :: set_id => cell_set_id
procedure :: universe => cell_universe
procedure :: fill => cell_fill
procedure :: material_size => cell_material_size
procedure :: material => cell_material
procedure :: sqrtkT_size => cell_sqrtkT_size
procedure :: sqrtkT => cell_sqrtkT
end type Cell
! array index of the root universe
integer(C_INT), bind(C) :: root_universe
integer(C_INT32_T), bind(C) :: n_cells ! # of cells
integer(C_INT32_T), bind(C) :: n_universes ! # of universes
type(Cell), allocatable, target :: cells(:)
! Dictionaries which map user IDs to indices in the global arrays
type(DictIntInt) :: cell_dict
contains
!===============================================================================
function cell_id(this) result(id)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: id
id = cell_id_c(this % ptr)
end function cell_id
subroutine cell_set_id(this, id)
class(Cell), intent(in) :: this
integer(C_INT32_T), intent(in) :: id
call cell_set_id_c(this % ptr, id)
end subroutine cell_set_id
function cell_universe(this) result(universe)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: universe
universe = cell_universe_c(this % ptr)
end function cell_universe
function cell_fill(this) result(fill)
class(Cell), intent(in) :: this
integer(C_INT32_T) :: fill
fill = cell_fill_c(this % ptr)
end function cell_fill
function cell_material_size(this) result(n)
class(Cell), intent(in) :: this
integer(C_INT) :: n
n = cell_material_size_c(this % ptr)
end function cell_material_size
function cell_material(this, i) result(mat)
class(Cell), intent(in) :: this
integer, intent(in) :: i
integer(C_INT32_T) :: mat
mat = cell_material_c(this % ptr, i)
end function cell_material
function cell_sqrtkT_size(this) result(n)
class(Cell), intent(in) :: this
integer :: n
n = cell_sqrtkT_size_c(this % ptr)
end function cell_sqrtkT_size
function cell_sqrtkT(this, i) result(sqrtkT)
class(Cell), intent(in) :: this
integer, intent(in) :: i
real(C_DOUBLE) :: sqrtkT
sqrtkT = cell_sqrtkT_c(this % ptr, i)
end function cell_sqrtkT
!===============================================================================
! FREE_MEMORY_GEOMETRY deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_geometry()
interface
subroutine free_memory_geometry_c() bind(C)
end subroutine free_memory_geometry_c
end interface
call free_memory_geometry_c()
n_cells = 0
n_universes = 0
if (allocated(cells)) deallocate(cells)
call cell_dict % clear()
end subroutine free_memory_geometry
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_cells(n, index_start, index_end) result(err) bind(C)
! Extend the cells array by n elements
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT32_T) :: i
integer(C_INT) :: err
type(Cell), allocatable :: temp(:) ! temporary cells array
if (n_cells == 0) then
! Allocate cells array
allocate(cells(n))
else
! Allocate cells array with increased size
allocate(temp(n_cells + n))
! Copy original cells to temporary array
temp(1:n_cells) = cells
! Move allocation from temporary array
call move_alloc(FROM=temp, TO=cells)
end if
! Return indices in cells array
if (present(index_start)) index_start = n_cells + 1
if (present(index_end)) index_end = n_cells + n
! Extend the C++ cells array and get pointers to the C++ objects
call extend_cells_c(n)
do i = n_cells - n, n_cells
cells(i) % ptr = cell_pointer(i - 1)
end do
err = 0
end function openmc_extend_cells
function openmc_get_cell_index(id, index) result(err) bind(C)
! Return the index in the cells array of a cell with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has(id)) then
index = cell_dict % get(id)
err = 0
else
err = E_INVALID_ID
call set_errmsg("No cell exists with ID=" // trim(to_str(id)) // ".")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory has not been allocated for cells.")
end if
end function openmc_get_cell_index
function openmc_cell_get_id(index, id) result(err) bind(C)
! Return the ID of a cell
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(cells)) then
id = cells(index) % id()
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_get_id
function openmc_cell_set_id(index, id) result(err) bind(C)
! Set the ID of a cell
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= n_cells) then
call cells(index) % set_id(id)
call cell_dict % set(id, index)
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in cells array is out of bounds.")
end if
end function openmc_cell_set_id
end module geometry_header

1
src/initialize.cpp
View file

@ -31,7 +31,6 @@
// data/functions from Fortran side
extern "C" void read_command_line();
extern "C" void read_geometry_xml();
extern "C" void read_materials_xml();
extern "C" void read_plots_xml();
extern "C" void read_tallies_xml();

312
src/input_xml.F90
View file

@ -2,24 +2,17 @@ module input_xml
use, intrinsic :: ISO_C_BINDING
use algorithm, only: find
use constants
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message, openmc_err_msg
use geometry_header
#ifdef DAGMC
use dagmc_header
#endif
use hdf5_interface
use material_header
use mesh_header
use message_passing
use mgxs_interface
use nuclide_header
use photon_header
use random_lcg, only: prn
use surface_header
use settings
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, &
@ -35,42 +28,11 @@ module input_xml
save
interface
subroutine count_cell_instances(univ_indx) bind(C)
import C_INT32_T
integer(C_INT32_T), intent(in), value :: univ_indx
end subroutine count_cell_instances
subroutine read_surfaces(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_surfaces
subroutine read_cells(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_cells
subroutine read_lattices(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_lattices
subroutine read_materials(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_materials
function find_root_universe() bind(C) result(root)
import C_INT32_T
integer(C_INT32_T) :: root
end function find_root_universe
function maximum_levels(univ) bind(C) result(n)
import C_INT32_T, C_INT
integer(C_INT32_T), intent(in), value :: univ
integer(C_INT) :: n
end function maximum_levels
subroutine read_plots(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
@ -86,268 +48,6 @@ module input_xml
contains
#ifdef DAGMC
!===============================================================================
! READ_GEOMETRY_DAGMC reads data from a DAGMC .h5m file, checking
! for material properties and surface boundary conditions
! some universe information is spoofed for now
!===============================================================================
subroutine read_geometry_dagmc()
integer :: i
integer :: univ_id
integer :: n_cells_in_univ
logical :: file_exists
character(MAX_LINE_LEN) :: filename
type(Cell), pointer :: c
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
! Check if dagmc.h5m exists
filename = trim(path_input) // "dagmc.h5m"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Geometry DAGMC file '" // trim(filename) // "' does not &
&exist!")
end if
call write_message("Reading DAGMC geometry...", 5)
call load_dagmc_geometry()
call allocate_surfaces()
call allocate_cells()
! setup universe data structs
do i = 1, n_cells
c => cells(i)
! additional metadata spoofing
univ_id = c % universe()
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
root_universe = find_root_universe()
end subroutine read_geometry_dagmc
#endif
!===============================================================================
! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_geometry_xml() bind(C)
integer :: i, n
integer :: univ_id
integer :: n_cells_in_univ
real(8) :: phi, theta, psi
logical :: file_exists
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
type(Cell), pointer :: c
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_cell
type(XMLNode), allocatable :: node_cell_list(:)
type(VectorInt) :: univ_ids ! List of all universe IDs
type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each
! universe contains
#ifdef DAGMC
if (dagmc) then
call read_geometry_dagmc()
return
end if
#endif
! Display output message
call write_message("Reading geometry XML file...", 5)
! Check if geometry.xml exists
filename = trim(path_input) // "geometry.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Geometry XML file '" // trim(filename) // "' does not &
&exist!")
end if
! Parse geometry.xml file
call doc % load_file(filename)
root = doc % document_element()
! ==========================================================================
! READ SURFACES FROM GEOMETRY.XML
! This variable is used to check whether at least one boundary condition was
! applied to a surface
boundary_exists = .false.
call read_surfaces(root % ptr)
! Allocate surfaces array
allocate(surfaces(surfaces_size()))
do i = 1, size(surfaces)
surfaces(i) % ptr = surface_pointer(i - 1);
if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
end do
! Check to make sure a boundary condition was applied to at least one
! surface
if (run_mode /= MODE_PLOTTING) then
if (.not. boundary_exists) then
call fatal_error("No boundary conditions were applied to any surfaces!")
end if
end if
! ==========================================================================
! READ CELLS FROM GEOMETRY.XML
call read_cells(root % ptr)
! Get pointer to list of XML <cell>
call get_node_list(root, "cell", node_cell_list)
! Get number of <cell> tags
n_cells = size(node_cell_list)
! Check for no cells
if (n_cells == 0) then
call fatal_error("No cells found in geometry.xml!")
end if
! Allocate cells array
allocate(cells(n_cells))
n_universes = 0
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer(i - 1)
! Get pointer to i-th cell node
node_cell = node_cell_list(i)
! Check to make sure 'id' hasn't been used
if (cell_dict % has(c % id())) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id()))
end if
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
! Rotations can only be applied to cells that are being filled with
! another universe
if (c % fill() == C_NONE) then
call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
&c % id())) // " because it is not filled with another universe")
end if
! Read number of rotation parameters
n = node_word_count(node_cell, "rotation")
if (n /= 3) then
call fatal_error("Incorrect number of rotation parameters on cell " &
// to_str(c % id()))
end if
! Copy rotation angles in x,y,z directions
allocate(c % rotation(3))
call get_node_array(node_cell, "rotation", c % rotation)
phi = -c % rotation(1) * PI/180.0_8
theta = -c % rotation(2) * PI/180.0_8
psi = -c % rotation(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
c % rotation_matrix = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
univ_id = c % universe()
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
! ==========================================================================
! READ LATTICES FROM GEOMETRY.XML
call read_lattices(root % ptr)
! ==========================================================================
! SETUP UNIVERSES
! Allocate universes, universe cell arrays, and assign base universe
root_universe = find_root_universe()
! Clear dictionary
call cells_in_univ_dict%clear()
! Close geometry XML file
call doc % clear()
end subroutine read_geometry_xml
subroutine allocate_surfaces()
integer :: i
! Allocate surfaces array
allocate(surfaces(surfaces_size()))
do i = 1, size(surfaces)
surfaces(i) % ptr = surface_pointer(i - 1);
end do
end subroutine allocate_surfaces
subroutine allocate_cells()
integer :: i
type(Cell), pointer :: c
! Allocate cells array
allocate(cells(n_cells))
do i = 1, n_cells
c => cells(i)
c % ptr = cell_pointer(i - 1)
! Check to make sure 'id' hasn't been used
if (cell_dict % has(c % id())) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id()))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)
end do
end subroutine allocate_cells
subroutine read_materials_xml() bind(C)
logical :: file_exists ! does materials.xml exist?
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
@ -355,11 +55,6 @@ contains
type(XMLNode) :: root
interface
function nuclides_size() bind(C) result(n)
import C_INT
integer(C_INT) :: n
end function
function elements_size() bind(C) result(n)
import C_INT
integer(C_INT) :: n
@ -395,9 +90,7 @@ contains
call read_materials(root % ptr)
! Set total number of nuclides and elements
n_nuclides = nuclides_size()
n_elements = elements_size()
allocate(nuclides(n_nuclides))
! Close materials XML file
call doc % clear()
@ -468,6 +161,11 @@ contains
import C_PTR
type(C_PTR) :: node_ptr
end subroutine
subroutine read_meshes(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine
end interface
! Check if tallies.xml exists

7
src/lattice.cpp
View file

@ -864,13 +864,12 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
// Non-method functions
//==============================================================================
extern "C" void
read_lattices(pugi::xml_node *node)
void read_lattices(pugi::xml_node node)
{
for (pugi::xml_node lat_node : node->children("lattice")) {
for (pugi::xml_node lat_node : node.children("lattice")) {
model::lattices.push_back(new RectLattice(lat_node));
}
for (pugi::xml_node lat_node : node->children("hex_lattice")) {
for (pugi::xml_node lat_node : node.children("hex_lattice")) {
model::lattices.push_back(new HexLattice(lat_node));
}

4
src/main.cpp
View file

@ -5,6 +5,7 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/particle_restart.h"
#include "openmc/settings.h"
@ -36,7 +37,8 @@ int main(int argc, char* argv[]) {
err = openmc_plot_geometry();
break;
case RUN_MODE_PARTICLE:
if (mpi::master) err = openmc_particle_restart();
if (mpi::master) run_particle_restart();
err = 0;
break;
case RUN_MODE_VOLUME:
err = openmc_calculate_volumes();

5
src/material.cpp
View file

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

8
src/material_header.F90
View file

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

20
src/math.F90
View file

@ -1,20 +0,0 @@
module math
use, intrinsic :: ISO_C_BINDING
implicit none
private
public :: t_percentile
interface
pure function t_percentile(p, df) bind(C) result(t)
use ISO_C_BINDING
implicit none
real(C_DOUBLE), value, intent(in) :: p
integer(C_INT), value, intent(in) :: df
real(C_DOUBLE) :: t
end function t_percentile
end interface
end module math

296
src/mesh_header.F90
View file

@ -1,296 +0,0 @@
module mesh_header
use, intrinsic :: ISO_C_BINDING
implicit none
!===============================================================================
! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
! congruent squares or cubes
!===============================================================================
type :: RegularMesh
type(C_PTR) :: ptr
contains
procedure :: id => regular_id
procedure :: volume_frac => regular_volume_frac
procedure :: n_dimension => regular_n_dimension
procedure :: dimension => regular_dimension
procedure :: lower_left => regular_lower_left
procedure :: upper_right => regular_upper_right
procedure :: width => regular_width
procedure :: get_bin => regular_get_bin
procedure :: get_indices => regular_get_indices
procedure :: get_bin_from_indices => regular_get_bin_from_indices
procedure :: get_indices_from_bin => regular_get_indices_from_bin
end type RegularMesh
interface
function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT) :: err
end function openmc_extend_meshes
function openmc_get_mesh_index(id, index) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
end function openmc_get_mesh_index
function openmc_mesh_get_id(index, id) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
end function openmc_mesh_get_id
function openmc_mesh_set_id(index, id) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
end function openmc_mesh_set_id
function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C)
import C_INT32_T, C_PTR, C_INT
integer(C_INT32_T), value, intent(in) :: index
type(C_PTR), intent(out) :: dims
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
end function openmc_mesh_get_dimension
function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
integer(C_INT), intent(in) :: dims(n)
integer(C_INT) :: err
end function openmc_mesh_set_dimension
function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C)
import C_INT32_T, C_PTR, C_INT
integer(C_INT32_T), value, intent(in) :: index
type(C_PTR), intent(out) :: ll
type(C_PTR), intent(out) :: ur
type(C_PTR), intent(out) :: width
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
end function openmc_mesh_get_params
function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C)
import C_INT32_T, C_INT, C_DOUBLE
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT), value, intent(in) :: n
real(C_DOUBLE), intent(in), optional :: ll(n)
real(C_DOUBLE), intent(in), optional :: ur(n)
real(C_DOUBLE), intent(in), optional :: width(n)
integer(C_INT) :: err
end function openmc_mesh_set_params
function mesh_id(ptr) result(id) bind(C)
import C_PTR, C_INT32_T
type(C_PTR), value :: ptr
integer(C_INT32_T) :: id
end function
function mesh_volume_frac(ptr) result(volume_frac) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE) :: volume_frac
end function
function mesh_n_dimension(ptr) result(n) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT) :: n
end function
function mesh_dimension(ptr, i) result(d) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: i
integer(C_INT) :: d
end function
function mesh_lower_left(ptr, i) result(ll) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: i
real(C_DOUBLE) :: ll
end function
function mesh_upper_right(ptr, i) result(ur) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: i
real(C_DOUBLE) :: ur
end function
function mesh_width(ptr, i) result(w) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: i
real(C_DOUBLE) :: w
end function
pure function mesh_get_bin(ptr, xyz) result(bin) bind(C)
import C_PTR, C_DOUBLE, C_INT
type(C_PTR), intent(in), value :: ptr
real(C_DOUBLE), intent(in) :: xyz(*)
integer(C_INT) :: bin
end function
pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in) :: ijk(*)
integer(C_INT) :: bin
end function
pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C)
import C_PTR, C_DOUBLE, C_INT, C_BOOL
type(C_PTR), intent(in), value :: ptr
real(C_DOUBLE), intent(in) :: xyz(*)
integer(C_INT), intent(out) :: ijk(*)
logical(C_BOOL), intent(out) :: in_mesh
end subroutine
pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in), value :: bin
integer(C_INT), intent(out) :: ijk(*)
end subroutine
function mesh_ptr(i) result(ptr) bind(C)
import C_INT, C_PTR
integer(C_INT), value :: i
type(C_PTR) :: ptr
end function
subroutine read_meshes(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine
function n_meshes() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
end interface
contains
function meshes(i) result(m)
integer, intent(in) :: i
type(RegularMesh) :: m
m % ptr = mesh_ptr(i)
end function
function regular_id(this) result(id)
class(RegularMesh), intent(in) :: this
integer(C_INT32_T) :: id
id = mesh_id(this % ptr)
end function
function regular_volume_frac(this) result(volume_frac)
class(RegularMesh), intent(in) :: this
real(C_DOUBLE) :: volume_frac
volume_frac = mesh_volume_frac(this % ptr)
end function
function regular_n_dimension(this) result(n)
class(RegularMesh), intent(in) :: this
integer(C_INT) :: n
n = mesh_n_dimension(this % ptr)
end function
function regular_dimension(this, i) result(d)
class(RegularMesh), intent(in) :: this
integer(C_INT), intent(in) :: i
integer(C_INT) :: d
d = mesh_dimension(this % ptr, i)
end function
function regular_lower_left(this, i) result(ll)
class(RegularMesh), intent(in) :: this
integer(C_INT), intent(in) :: i
real(C_DOUBLE) :: ll
ll = mesh_lower_left(this % ptr, i)
end function
function regular_upper_right(this, i) result(ur)
class(RegularMesh), intent(in) :: this
integer(C_INT), intent(in) :: i
real(C_DOUBLE) :: ur
ur = mesh_upper_right(this % ptr, i)
end function
function regular_width(this, i) result(w)
class(RegularMesh), intent(in) :: this
integer(C_INT), intent(in) :: i
real(C_DOUBLE) :: w
w = mesh_width(this % ptr, i)
end function
!===============================================================================
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
pure subroutine regular_get_bin(this, xyz, bin)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz(:) ! coordinates
integer, intent(out) :: bin ! tally bin
bin = mesh_get_bin(this % ptr, xyz)
end subroutine
!===============================================================================
! GET_MESH_INDICES determines the indices of a particle in a structured mesh
!===============================================================================
pure subroutine regular_get_indices(this, xyz, ijk, in_mesh)
class(RegularMesh), intent(in) :: this
real(8), intent(in) :: xyz(:) ! coordinates to check
integer, intent(out) :: ijk(:) ! indices in mesh
logical, intent(out) :: in_mesh ! were given coords in mesh?
logical(C_BOOL) :: in_mesh_
call mesh_get_indices(this % ptr, xyz, ijk, in_mesh_)
in_mesh = in_mesh_
end subroutine regular_get_indices
!===============================================================================
! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number
! for use in a TallyObject results array
!===============================================================================
pure function regular_get_bin_from_indices(this, ijk) result(bin)
class(RegularMesh), intent(in) :: this
integer, intent(in) :: ijk(:)
integer :: bin
bin = mesh_get_bin_from_indices(this % ptr, ijk)
end function regular_get_bin_from_indices
!===============================================================================
! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to
! (i), (i,j), or (i,j,k) indices
!===============================================================================
pure subroutine regular_get_indices_from_bin(this, bin, ijk)
class(RegularMesh), intent(in) :: this
integer, intent(in) :: bin
integer, intent(out) :: ijk(:)
call mesh_get_indices_from_bin(this % ptr, bin, ijk)
end subroutine regular_get_indices_from_bin
end module mesh_header

82
src/nuclide.cpp
View file

@ -892,18 +892,11 @@ void check_data_version(hid_t file_id)
// C API
//==============================================================================
extern "C" void extend_nuclides();
extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr,
const double* temps, int n, int n_nuclide);
extern "C" int openmc_load_nuclide(const char* name)
{
if (data::nuclide_map.find(name) == data::nuclide_map.end()) {
const auto& it = data::library_map.find({Library::Type::neutron, name});
if (it != data::library_map.end()) {
// Extend nuclides array on Fortran side
extend_nuclides();
// Get filename for library containing nuclide
int idx = it->second;
std::string& filename = data::libraries[idx].path_;
@ -920,10 +913,6 @@ extern "C" int openmc_load_nuclide(const char* name)
data::nuclides.push_back(std::make_unique<Nuclide>(
group, temperature, i_nuclide));
// Read from Fortran too
nuclide_from_hdf5(group, data::nuclides.back().get(),
&temperature.front(), temperature.size(), i_nuclide + 1);
close_group(group);
file_close(file_id);
@ -943,6 +932,30 @@ extern "C" int openmc_load_nuclide(const char* name)
return 0;
}
extern "C" int
openmc_get_nuclide_index(const char* name, int* index)
{
auto it = data::nuclide_map.find(name);
if (it == data::nuclide_map.end()) {
set_errmsg("No nuclide named '" + std::string{name} + "' has been loaded.");
return OPENMC_E_DATA;
}
*index = it->second;
return 0;
}
extern "C" int
openmc_nuclide_name(int index, const char** name)
{
if (index >= 0 && index < data::nuclides.size()) {
*name = data::nuclides[index]->name_.data();
return 0;
} else {
set_errmsg("Index in nuclides vector is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -956,47 +969,6 @@ set_particle_energy_bounds(int particle, double E_min, double E_max)
extern "C" int nuclides_size() { return data::nuclide_map.size(); }
extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); }
extern "C" double nuclide_awr(int i_nuc) { return data::nuclides[i_nuc - 1]->awr_; }
extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx)
{
return nuc->reactions_[i_rx-1].get();
}
extern "C" void nuclide_calculate_xs_c(Nuclide* nuc, int i_sab, double E,
int i_log_union, double sqrtkT, double sab_frac)
{
nuc->calculate_xs(i_sab, E, i_log_union, sqrtkT, sab_frac);
}
extern "C" void nuclide_calculate_elastic_xs_c(Nuclide* nuc)
{
nuc->calculate_elastic_xs();
}
extern "C" double nuclide_nu_c(Nuclide* nuc, double E, int emission_mode, int group)
{
return nuc->nu(E, static_cast<Nuclide::EmissionMode>(emission_mode - 1), group);
}
extern "C" double nuclide_fission_q_prompt(Nuclide* nuc, double E)
{
return nuc->fission_q_prompt_ ? (*nuc->fission_q_prompt_)(E) : 0.0;
}
extern "C" double nuclide_fission_q_recov(Nuclide* nuc, double E)
{
return nuc->fission_q_recov_ ? (*nuc->fission_q_recov_)(E) : 0.0;
}
extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT,
double* sig_s, double* sig_a, double* sig_f)
{
std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT);
}
extern "C" bool multipole_in_range(const Nuclide* nuc, double E)
{
return nuc->multipole_ && E >= nuc->multipole_->E_min_&&
@ -1009,12 +981,6 @@ extern "C" void nuclides_clear()
data::nuclide_map.clear();
}
extern "C" int nuclide_map_get(const char* name)
{
auto it = data::nuclide_map.find(name);
return it == data::nuclide_map.end() ? -1 : it->second + 1;
}
extern "C" NuclideMicroXS* micro_xs_ptr();
extern "C" ElementMicroXS* micro_photon_xs_ptr();

615
src/nuclide_header.F90
View file

@ -3,82 +3,10 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use algorithm, only: sort, find
use constants
use endf, only: is_fission, is_disappearance
use error
use hdf5_interface
use message_passing
use reaction_header, only: Reaction
use settings
use stl_vector, only: VectorInt, VectorReal
use string
implicit none
!===============================================================================
! Nuclide contains the base nuclidic data for a nuclide described as needed
! for continuous-energy neutron transport.
!===============================================================================
type EnergyGrid
real(8), allocatable :: energy(:) ! energy values corresponding to xs
end type EnergyGrid
! Positions for first dimension of Nuclide % xs
integer, parameter :: &
XS_TOTAL = 1, &
XS_ABSORPTION = 2, &
XS_FISSION = 3, &
XS_NU_FISSION = 4, &
XS_PHOTON_PROD = 5
! The array within SumXS is of shape (5, n_energy) where the first dimension
! corresponds to the following values: 1) total, 2) absorption (MT > 100), 3)
! fission, 4) neutron production, 5) photon production
type SumXS
real(8), allocatable :: value(:,:)
end type SumXS
type :: Nuclide
! Nuclide meta-data
character(20) :: name ! name of nuclide, e.g. U235
integer :: Z ! atomic number
integer :: A ! mass number
integer :: metastable ! metastable state
real(8) :: awr ! Atomic Weight Ratio
integer :: i_nuclide ! The nuclides index in the nuclides array
real(8), allocatable :: kTs(:) ! temperature in eV (k*T)
! Fission information
logical :: fissionable = .false. ! nuclide is fissionable?
! Energy grid for each temperature
type(EnergyGrid), allocatable :: grid(:)
! Microscopic cross sections
type(SumXS), allocatable :: xs(:)
! Fission information
integer, allocatable :: index_fission(:) ! indices in reactions
! Reactions
type(Reaction), allocatable :: reactions(:)
! Array that maps MT values to index in reactions; used at tally-time. Note
! that ENDF-102 does not assign any MT values above 891.
integer :: reaction_index(891)
type(C_PTR) :: ptr
contains
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
procedure :: calculate_xs => nuclide_calculate_xs
procedure :: calculate_elastic_xs => nuclide_calculate_elastic_xs
end type Nuclide
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a particular
! nuclide at the current energy
@ -103,7 +31,7 @@ module nuclide_header
! Cross sections for depletion reactions (note that these are not stored in
! macroscopic cache)
real(C_DOUBLE) :: reaction(size(DEPLETION_RX))
real(C_DOUBLE) :: reaction(6)
! Indicies and factors needed to compute cross sections from the data tables
integer(C_INT) :: index_grid ! Index on nuclide energy grid
@ -140,487 +68,25 @@ module nuclide_header
real(C_DOUBLE) :: pair_production ! macroscopic pair production xs
end type MaterialMacroXS
! Nuclear data for each nuclide
type(Nuclide), allocatable, target :: nuclides(:)
integer(C_INT), bind(C) :: n_nuclides
! Cross section caches
type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide
type(MaterialMacroXS), bind(C) :: material_xs ! Cache for current material
!$omp threadprivate(micro_xs, material_xs)
interface
function library_present_c(type, name) result(b) bind(C, name='library_present')
import C_INT, C_CHAR, C_BOOL
integer(C_INT), value :: type
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: b
end function
function library_path_c(type, name) result(path) bind(C, name='library_path')
import C_INT, C_CHAR, C_PTR
integer(C_INT), value :: type
character(kind=C_CHAR), intent(in) :: name(*)
type(C_PTR) :: path
end function
subroutine multipole_deriv_eval(ptr, E, sqrtkT, sig_s, sig_a, sig_f) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE), value :: sqrtkT
real(C_DOUBLE), intent(out) :: sig_s
real(C_DOUBLE), intent(out) :: sig_a
real(C_DOUBLE), intent(out) :: sig_f
end subroutine
function multipole_in_range(ptr, E) result(b) bind(C)
import C_PTR, C_DOUBLE, C_BOOL
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
logical(C_BOOL) :: b
end function
function nuclide_awr(i_nuc) result(awr) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: i_nuc
real(C_DOUBLE) :: awr
end function
function nuclide_fission_q_prompt(ptr, E) result(q) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE) :: q
end function
function nuclide_fission_q_recov(ptr, E) result(q) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE) :: q
end function
function nuclide_map_get(name) result(idx) bind(C)
import C_CHAR, C_INT
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT) :: idx
function nuclides_size() bind(C) result(n)
import C_INT
integer(C_INT) :: n
end function
end interface
contains
function library_path(type, name) result(path)
integer, intent(in) :: type
character(len=*), intent(in) :: name
character(MAX_FILE_LEN) :: path
type(C_PTR) :: ptr
character(kind=C_CHAR), pointer :: string(:)
ptr = library_path_c(type, to_c_string(name))
call c_f_pointer(ptr, string, [255])
path = to_f_string(string)
end function
function library_present(type, name) result(b)
integer, intent(in) :: type
character(len=*), intent(in) :: name
logical :: b
b = library_present_c(type, to_c_string(name))
end function
function micro_xs_ptr() result(ptr) bind(C)
type(C_PTR) :: ptr
ptr = C_LOC(micro_xs(1))
end function
subroutine nuclide_from_hdf5(group_id, ptr, temps, n, i_nuclide) bind(C)
integer(HID_T), value :: group_id
type(C_PTR), value :: ptr
type(C_PTR), value :: temps
integer(C_INT), value :: n
integer(C_INT), value :: i_nuclide
real(C_DOUBLE), pointer :: temperature(:) ! list of desired temperatures
integer :: i
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: energy_group, energy_dset
integer(HID_T) :: kT_group
integer(HID_T) :: rxs_group
integer(HID_T) :: rx_group
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN), allocatable :: dset_names(:)
character(MAX_WORD_LEN), allocatable :: grp_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
! Get array passed
call c_f_pointer(temps, temperature, [n])
associate (this => nuclides(i_nuclide))
this % ptr = ptr
! Get name of nuclide from group
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))
call read_attribute(this % Z, group_id, 'Z')
call read_attribute(this % A, group_id, 'A')
call read_attribute(this % metastable, group_id, 'metastable')
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
kT_group = open_group(group_id, 'kTs')
! Determine temperatures available
call get_datasets(kT_group, dset_names)
allocate(temps_available(size(dset_names)))
do i = 1, size(dset_names)
! Read temperature value
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
temps_available(i) = temps_available(i) / K_BOLTZMANN
end do
call sort(temps_available)
! If only one temperature is available, revert to nearest temperature
if (size(temps_available) == 1 .and. temperature_method == TEMPERATURE_INTERPOLATION) then
if (master) then
call warning("Cross sections for " // trim(this % name) // " are only &
&available at one temperature. Reverting to nearest temperature &
&method.")
end if
temperature_method = TEMPERATURE_NEAREST
end if
! Determine actual temperatures to read -- start by checking whether a
! temperature range was given, in which case all temperatures in the range
! are loaded irrespective of what temperatures actually appear in the model
if (temperature_range(2) > ZERO) then
do i = 1, size(temps_available)
temp_actual = temps_available(i)
if (temperature_range(1) <= temp_actual .and. temp_actual <= temperature_range(2)) then
call temps_to_read % push_back(nint(temp_actual))
end if
end do
end if
select case (temperature_method)
case (TEMPERATURE_NEAREST)
! Find nearest temperatures
do i = 1, n
temp_desired = temperature(i)
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
temp_actual = temps_available(i_closest)
if (abs(temp_actual - temp_desired) < temperature_tolerance) then
if (find(temps_to_read, nint(temp_actual)) == -1) then
call temps_to_read % push_back(nint(temp_actual))
! Write warning for resonance scattering data if 0K is not available
if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0 &
.and. master) then
call warning(trim(this % name) // " does not contain 0K data &
&needed for resonance scattering options selected. Using &
&data at " // trim(to_str(temp_actual)) &
// " K instead.")
end if
end if
else
call fatal_error("Nuclear data library does not contain cross &
&sections for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
end if
end do
case (TEMPERATURE_INTERPOLATION)
! If temperature interpolation or multipole is selected, get a list of
! bounding temperatures for each actual temperature present in the model
TEMP_LOOP: do i = 1, n
temp_desired = temperature(i)
do j = 1, size(temps_available) - 1
if (temps_available(j) <= temp_desired .and. &
temp_desired < temps_available(j + 1)) then
if (find(temps_to_read, nint(temps_available(j))) == -1) then
call temps_to_read % push_back(nint(temps_available(j)))
end if
if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then
call temps_to_read % push_back(nint(temps_available(j + 1)))
end if
cycle TEMP_LOOP
end if
end do
call fatal_error("Nuclear data library does not contain cross sections &
&for " // trim(this % name) // " at temperatures that bound " // &
trim(to_str(nint(temp_desired))) // " K.")
end do TEMP_LOOP
end select
! Sort temperatures to read
call sort(temps_to_read)
n_temperature = temps_to_read % size()
allocate(this % kTs(n_temperature))
allocate(this % grid(n_temperature))
! Get kT values
do i = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
! Read exact temperature value
call read_dataset(this % kTs(i), kT_group, trim(temp_str))
end do
call close_group(kT_group)
! Read energy grid
energy_group = open_group(group_id, 'energy')
do i = 1, n_temperature
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
energy_dset = open_dataset(energy_group, temp_str)
call get_shape(energy_dset, dims)
allocate(this % grid(i) % energy(int(dims(1), 4)))
call read_dataset(this % grid(i) % energy, energy_dset)
call close_dataset(energy_dset)
end do
call close_group(energy_group)
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
call get_groups(rxs_group, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "reaction_")) then
call MTs % push_back(int(str_to_int(grp_names(j)(10:12))))
end if
end do
! Read reactions
allocate(this % reactions(MTs % size()))
do i = 1, size(this % reactions)
rx_group = open_group(rxs_group, 'reaction_' // trim(&
zero_padded(MTs % data(i), 3)))
! Set pointer for each reaction
call this % reactions(i) % init(this % ptr, i)
call close_group(rx_group)
end do
call close_group(rxs_group)
! Create derived cross section data
call this % create_derived()
! Finalize with the nuclide index
this % i_nuclide = i_nuclide
end associate
end subroutine nuclide_from_hdf5
subroutine nuclide_create_derived(this)
class(Nuclide), intent(inout) :: this
integer :: i, j, k, l
integer :: t
integer :: n
integer :: n_grid
integer :: i_fission
integer :: n_temperature
type(VectorInt) :: MTs
n_temperature = size(this % kTs)
allocate(this % xs(n_temperature))
this % reaction_index(:) = 0
do i = 1, n_temperature
! Allocate and initialize derived cross sections
n_grid = size(this % grid(i) % energy)
allocate(this % xs(i) % value(5,n_grid))
this % xs(i) % value(:,:) = ZERO
end do
i_fission = 0
do i = 1, size(this % reactions)
call MTs % push_back(this % reactions(i) % MT)
this % reaction_index(this % reactions(i) % MT) = i
associate (rx => this % reactions(i))
do t = 1, n_temperature
j = rx % xs_threshold(t)
n = rx % xs_size(t)
! Calculate photon production cross section
do k = 1, rx % products_size()
if (rx % product_particle(k) == PHOTON) then
do l = 1, n
this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = &
this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + &
rx % xs(t, l) * rx % product_yield(k, &
this % grid(t) % energy(l+j-1))
end do
end if
end do
! Skip gas production cross sections (MT=200+), etc.
if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle
! Skip any reaction that has been marked as redundant
if (rx % redundant) cycle
! Add contribution to total cross section
do k = j, j + n - 1
this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % &
value(XS_TOTAL,k) + rx % xs(t, k - j + 1)
end do
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
do k = j, j + n - 1
this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % &
value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1)
end do
end if
! Information about fission reactions
if (t == 1) then
if (rx % MT == N_FISSION) then
allocate(this % index_fission(1))
elseif (rx % MT == N_F) then
allocate(this % index_fission(PARTIAL_FISSION_MAX))
end if
end if
! Add contribution to fission cross section
if (is_fission(rx % MT)) then
this % fissionable = .true.
do k = j, j + n - 1
this % xs(t) % value(XS_FISSION,k) = this % xs(t) % &
value(XS_FISSION,k) + rx % xs(t, k - j + 1)
! Also need to add fission cross sections to absorption
this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % &
value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1)
end do
! Keep track of this reaction for easy searching later
if (t == 1) then
i_fission = i_fission + 1
this % index_fission(i_fission) = i
end if
end if ! fission
end do ! temperature
end associate ! rx
end do ! reactions
! Calculate nu-fission cross section
do t = 1, n_temperature
if (this % fissionable) then
n_grid = size(this % grid(t) % energy)
do i = 1, n_grid
this % xs(t) % value(XS_NU_FISSION,i) = &
this % nu(this % grid(t) % energy(i), EMISSION_TOTAL) * &
this % xs(t) % value(XS_FISSION,i)
end do
end if
end do
end subroutine nuclide_create_derived
!===============================================================================
! NUCLIDE_NU is an interface to the number of fission neutrons produced
!===============================================================================
pure function nuclide_nu(this, E, emission_mode, group) result(nu)
class(Nuclide), intent(in) :: this
real(8), intent(in) :: E
integer, intent(in) :: emission_mode
integer, optional, intent(in) :: group
real(8) :: nu
interface
pure function nuclide_nu_c(ptr, E, emission_mode, group) result(nu) bind(C)
import C_PTR, C_DOUBLE, C_INT
type(C_PTR), value, intent(in) :: ptr
real(C_DOUBLE), value, intent(in) :: E
integer(C_INT), value, intent(in) :: emission_mode
integer(C_INT), value, intent(in) :: group
real(C_DOUBLE) :: nu
end function
end interface
if (present(group)) then
nu = nuclide_nu_c(this % ptr, E, emission_mode, group)
else
nu = nuclide_nu_c(this % ptr, E, emission_mode, 0)
end if
end function nuclide_nu
subroutine nuclide_init_grid(this)
class(nuclide), intent(inout) :: this
interface
subroutine nuclide_init_grid_c(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call nuclide_init_grid_c(this % ptr)
end subroutine
!===============================================================================
! NUCLIDE_CALCULATE_XS determines microscopic cross sections for the nuclide
! at the energy of the given particle
!===============================================================================
subroutine nuclide_calculate_xs(this, i_sab, E, i_log_union, sqrtkT, &
sab_frac)
class(Nuclide), intent(in) :: this ! Nuclide object
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! square root of kT, material dependent
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
interface
subroutine nuclide_calculate_xs_c(ptr, i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: i_sab
real(C_DOUBLE), value :: E
integer(C_INT), value :: i_log_union
real(C_DOUBLE), value :: sqrtkT
real(C_DOUBLE), value :: sab_frac
end subroutine
end interface
call nuclide_calculate_xs_c(this % ptr, i_sab - 1, E, i_log_union, sqrtkT, sab_frac)
end subroutine nuclide_calculate_xs
!===============================================================================
! NUCLIDE_CALCULATE_ELASTIC_XS precalculates the free atom elastic scattering
! cross section. Normally it is not needed until a collision actually occurs in
! a material. However, in the thermal and unresolved resonance regions, we have
! to calculate it early to adjust the total cross section correctly.
!===============================================================================
subroutine nuclide_calculate_elastic_xs(this)
class(Nuclide), intent(in) :: this
interface
subroutine nuclide_calculate_elastic_xs_c(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call nuclide_calculate_elastic_xs_c(this % ptr)
end subroutine nuclide_calculate_elastic_xs
!===============================================================================
! FREE_MEMORY_NUCLIDE deallocates global arrays defined in this module
@ -635,80 +101,9 @@ contains
end subroutine
end interface
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) then
! First call the clear routines
deallocate(nuclides)
call nuclides_clear()
end if
n_nuclides = 0
call nuclides_clear()
call library_clear()
end subroutine free_memory_nuclide
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_get_nuclide_index(name, index) result(err) bind(C)
! Return the index in the nuclides array of a nuclide with a given name
character(kind=C_CHAR), intent(in) :: name(*)
integer(C_INT), intent(out) :: index
integer(C_INT) :: err
character(:), allocatable :: name_
! Copy array of C_CHARs to normal Fortran string
name_ = to_f_string(name)
err = 0
if (allocated(nuclides)) then
index = nuclide_map_get(name)
if (index == -1) then
err = E_DATA
call set_errmsg("No nuclide named '" // trim(name_) // &
"' has been loaded.")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory for nuclides has not been allocated.")
end if
end function openmc_get_nuclide_index
function openmc_nuclide_name(index, name) result(err) bind(C)
! Return the name of a nuclide with a given index
integer(C_INT), value, intent(in) :: index
type(c_ptr), intent(out) :: name
integer(C_INT) :: err
character(C_CHAR), pointer :: name_
err = E_UNASSIGNED
if (allocated(nuclides)) then
if (index >= 1 .and. index <= size(nuclides)) then
name_ => nuclides(index) % name(1:1)
name = C_LOC(name_)
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in nuclides array is out of bounds.")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory for nuclides has not been allocated yet.")
end if
end function openmc_nuclide_name
subroutine extend_nuclides() bind(C)
type(Nuclide), allocatable :: new_nuclides(:)
! allocate extra space in nuclides array
allocate(new_nuclides(n_nuclides + 1))
new_nuclides(1:n_nuclides) = nuclides(:)
call move_alloc(FROM=new_nuclides, TO=nuclides)
n_nuclides = n_nuclides + 1
end subroutine
end module nuclide_header

4
src/output.cpp
View file

@ -65,7 +65,7 @@ void title()
" ############### %%%%%%%%%%%%%%%%\n" <<
" ############ %%%%%%%%%%%%%%%\n" <<
" ######## %%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%\n";
" %%%%%%%%%%%\n\n";
// Write version information
std::cout <<
@ -88,7 +88,7 @@ void title()
#ifdef _OPENMP
// Write number of OpenMP threads
std::cout << " OpenMC Threads | " << omp_get_max_threads() << '\n';
std::cout << " OpenMP Threads | " << omp_get_max_threads() << '\n';
#endif
std::cout << '\n';
}

512
src/particle.cpp
View file

@ -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,503 @@ 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
simulation::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
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;
} 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) {
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
auto surf_p = dynamic_cast<PeriodicSurface*>(surf);
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) {
auto cellp = dynamic_cast<DAGCell*>(model::cells[last_cell[0]]);
// TODO: off-by-one
auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface) - 1]);
int32_t i_cell = next_cell(cellp, surfp) - 1;
// save material and temp
last_material = material;
last_sqrtkT = sqrtkT;
// set new cell value
coord[0].cell = i_cell;
cell_instance = 0;
// TODO: off-by-one
int mat = model::cells[i_cell]->material_[0];
material = (mat == MATERIAL_VOID) ? mat : mat + 1;
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)
{

142
src/particle_restart.F90
View file

@ -1,142 +0,0 @@
module particle_restart
use, intrinsic :: ISO_FORTRAN_ENV
use bank_header, only: Bank
use constants
use error, only: write_message
use hdf5_interface, only: file_open, file_close, read_dataset, HID_T
use mgxs_interface, only: energy_bin_avg
use nuclide_header, only: micro_xs, n_nuclides
use particle_header
use photon_header, only: micro_photon_xs, n_elements
use random_lcg, only: set_particle_seed
use settings
use simulation_header
use tally_header, only: n_tallies
use tracking, only: transport
implicit none
private
public :: openmc_particle_restart
contains
!===============================================================================
! OPENMC_PARTICLE_RESTART is the main routine that runs the particle restart
!===============================================================================
function openmc_particle_restart() result(err) bind(C)
integer(C_INT) :: err
integer(8) :: particle_seed
integer :: previous_run_mode
type(Particle) :: p
interface
subroutine set_micro_xs() bind(C)
end subroutine
subroutine print_particle(p) bind(C)
import Particle
type(Particle), intent(in) :: p
end subroutine
end interface
err = 0
! Set verbosity high
verbosity = 10
!$omp parallel
allocate(micro_xs(n_nuclides))
allocate(micro_photon_xs(n_elements))
!$omp end parallel
call set_micro_xs()
! Initialize the particle to be tracked
call particle_initialize(p)
! Read in the restart information
call read_particle_restart(p, previous_run_mode)
! Set all tallies to 0 for now (just tracking errors)
n_tallies = 0
! Compute random number seed
select case (previous_run_mode)
case (MODE_EIGENVALUE)
particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id
case (MODE_FIXEDSOURCE)
particle_seed = p % id
end select
call set_particle_seed(particle_seed)
! Transport neutron
call transport(p)
! Write output if particle made it
call print_particle(p)
deallocate(micro_xs)
end function openmc_particle_restart
!===============================================================================
! READ_PARTICLE_RESTART reads the particle restart file
!===============================================================================
subroutine read_particle_restart(p, previous_run_mode)
type(Particle), intent(inout) :: p
integer, intent(inout) :: previous_run_mode
integer(HID_T) :: file_id
character(MAX_WORD_LEN) :: tempstr
! Write meessage
call write_message("Loading particle restart file " &
// trim(path_particle_restart) // "...", 5)
! Open file
file_id = file_open(path_particle_restart, 'r')
! Read data from file
call read_dataset(current_batch, file_id, 'current_batch')
call read_dataset(gen_per_batch, file_id, 'generations_per_batch')
call read_dataset(current_gen, file_id, 'current_generation')
call read_dataset(n_particles, file_id, 'n_particles')
call read_dataset(tempstr, file_id, 'run_mode')
select case (tempstr)
case ('eigenvalue')
previous_run_mode = MODE_EIGENVALUE
case ('fixed source')
previous_run_mode = MODE_FIXEDSOURCE
end select
call read_dataset(p % id, file_id, 'id')
call read_dataset(p % type, file_id, 'type')
call read_dataset(p % wgt, file_id, 'weight')
call read_dataset(p % E, file_id, 'energy')
call read_dataset(p % coord(1) % xyz, file_id, 'xyz')
call read_dataset(p % coord(1) % uvw, file_id, 'uvw')
! Set energy group and average energy in multi-group mode
if (.not. run_CE) then
p % g = int(p % E)
p % E = energy_bin_avg(p % g)
end if
! Set particle last attributes
p % last_wgt = p % wgt
p % last_xyz_current = p % coord(1)%xyz
p % last_xyz = p % coord(1)%xyz
p % last_uvw = p % coord(1)%uvw
p % last_E = p % E
p % last_g = p % g
! Close hdf5 file
call file_close(file_id)
end subroutine read_particle_restart
end module particle_restart

109
src/particle_restart.cpp Normal file
View file

@ -0,0 +1,109 @@
#include "openmc/particle_restart.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
#include <algorithm> // for copy
#include <array>
#include <string>
namespace openmc {
void read_particle_restart(Particle& p, int& previous_run_mode)
{
// Write meessage
write_message("Loading particle restart file " +
settings::path_particle_restart + "...", 5);
// Open file
hid_t file_id = file_open(settings::path_particle_restart, 'r');
// Read data from file
read_dataset(file_id, "current_batch", simulation::current_batch);
read_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
read_dataset(file_id, "current_generation", simulation::current_gen);
read_dataset(file_id, "n_particles", settings::n_particles);
std::string mode;
read_dataset(file_id, "run_mode", mode);
if (mode == "eigenvalue") {
previous_run_mode = RUN_MODE_EIGENVALUE;
} else if (mode == "fixed source") {
previous_run_mode = RUN_MODE_FIXEDSOURCE;
}
read_dataset(file_id, "id", p.id);
read_dataset(file_id, "type", p.type);
read_dataset(file_id, "weight", p.wgt);
read_dataset(file_id, "energy", p.E);
std::array<double, 3> x;
read_dataset(file_id, "xyz", x);
std::copy(x.data(), x.data() + 3, p.coord[0].xyz);
read_dataset(file_id, "uvw", x);
std::copy(x.data(), x.data() + 3, p.coord[0].uvw);
// Set energy group and average energy in multi-group mode
if (!settings::run_CE) {
p.g = p.E;
p.E = data::energy_bin_avg[p.g - 1];
}
// Set particle last attributes
p.last_wgt = p.wgt;
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current);
std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz);
std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw);
p.last_E = p.E;
p.last_g = p.g;
// Close hdf5 file
file_close(file_id);
}
void run_particle_restart()
{
// Set verbosity high
settings::verbosity = 10;
simulation_init_f();
set_micro_xs();
// Initialize the particle to be tracked
Particle p;
p.initialize();
// Read in the restart information
int previous_run_mode;
read_particle_restart(p, previous_run_mode);
// Set all tallies to 0 for now (just tracking errors)
model::tallies.clear();
// Compute random number seed
int64_t particle_seed;
switch (previous_run_mode) {
case RUN_MODE_EIGENVALUE:
particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id;
break;
case RUN_MODE_FIXEDSOURCE:
particle_seed = p.id;
break;
}
set_particle_seed(particle_seed);
// Transport neutron
p.transport();
// Write output if particle made it
print_particle(&p);
simulation_finalize_f();
}
} // namespace openmc

270
src/reaction_header.F90
View file

@ -1,270 +0,0 @@
module reaction_header
use, intrinsic :: ISO_C_BINDING
use hdf5_interface
implicit none
private
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type, public :: Reaction
type(C_PTR) :: ptr
integer(C_INT) :: MT ! ENDF MT value
real(C_DOUBLE) :: Q_value ! Reaction Q value
logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass?
logical(C_BOOL) :: redundant ! redundant reaction?
contains
procedure :: init
procedure :: mt_
procedure :: q_value_
procedure :: scatter_in_cm_
procedure :: redundant_
procedure :: product_decay_rate
procedure :: product_emission_mode
procedure :: product_particle
procedure :: product_sample
procedure :: product_yield
procedure :: products_size
procedure :: sample_elastic_mu
procedure :: xs
procedure :: xs_size
procedure :: xs_threshold
end type Reaction
interface
function nuclide_reaction(nuc_ptr, i_rx) result(ptr) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: nuc_ptr
integer(C_INT), value :: i_rx
type(C_PTR) :: ptr
end function
function reaction_mt(ptr) result(mt) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT) :: mt
end function
function reaction_q_value(ptr) result(q_value) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE) :: q_value
end function
function reaction_scatter_in_cm(ptr) result(b) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: ptr
logical(C_BOOL) :: b
end function
function reaction_redundant(ptr) result(b) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: ptr
logical(C_BOOL) :: b
end function
pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in), value :: product
real(C_DOUBLE) :: rate
end function
pure function reaction_product_emission_mode(ptr, product) result(m) bind(C)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in), value :: product
integer(C_INT) :: m
end function
pure function reaction_product_particle(ptr, product) result(particle) bind(C)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in), value :: product
integer(C_INT) :: particle
end function
subroutine reaction_product_sample(ptr, product, E_in, E_out, mu) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: product
real(C_DOUBLE), value :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
end subroutine
pure function reaction_product_yield(ptr, product, E) result(val) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), intent(in), value :: ptr
integer(C_INT), intent(in), value :: product
real(C_DOUBLE), intent(in), value :: E
real(C_DOUBLE) :: val
end function
pure function reaction_products_size(ptr) result(sz) bind(C)
import C_PTR, C_INT
type(C_PTR), intent(in), value :: ptr
integer(C_INT) :: sz
end function
function reaction_sample_elastic_mu(ptr, E) result(mu) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE) :: mu
end function
function reaction_xs(ptr, temperature, energy) result(xs) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT), value :: energy
real(C_DOUBLE) :: xs
end function
function reaction_xs_size(ptr, temperature) result(sz) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT) :: sz
end function
function reaction_xs_threshold(ptr, temperature) result(threshold) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT) :: threshold
end function
end interface
contains
subroutine init(this, nuc_ptr, i_rx)
class(Reaction), intent(inout) :: this
type(C_PTR), intent(in) :: nuc_ptr
integer(C_INT), intent(in) :: i_rx
this % ptr = nuclide_reaction(nuc_ptr, i_rx)
this % MT = reaction_mt(this % ptr)
this % Q_value = reaction_q_value(this % ptr)
this % scatter_in_cm = reaction_scatter_in_cm(this % ptr)
this % redundant = reaction_redundant(this % ptr)
end subroutine
function mt_(this) result(mt)
class(Reaction), intent(in) :: this
integer(C_INT) :: MT
mt = reaction_mt(this % ptr)
end function
function q_value_(this) result(q_value)
class(Reaction), intent(in) :: this
real(C_DOUBLE) :: q_value
q_value = reaction_q_value(this % ptr)
end function
function scatter_in_cm_(this) result(cm)
class (Reaction), intent(in) :: this
logical(C_BOOL) :: cm
cm = reaction_scatter_in_cm(this % ptr)
end function
function redundant_(this) result(redundant)
class (Reaction), intent(in) :: this
logical(C_BOOL) :: redundant
redundant = reaction_redundant(this % ptr)
end function
pure function product_decay_rate(this, product) result(rate)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE) :: rate
rate = reaction_product_decay_rate(this % ptr, product)
end function
pure function product_emission_mode(this, product) result(m)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
integer(C_INT) :: m
m = reaction_product_emission_mode(this % ptr, product)
end function
pure function product_particle(this, product) result(p)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
integer(C_INT) :: p
p = reaction_product_particle(this % ptr, product)
end function
subroutine product_sample(this, product, E_in, E_out, mu)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE), intent(in) :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
call reaction_product_sample(this % ptr, product, E_in, E_out, mu)
end subroutine
pure function product_yield(this, product, E) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE), intent(in) :: E
real(C_DOUBLE) :: val
val = reaction_product_yield(this % ptr, product, E)
end function
pure function products_size(this) result(sz)
class(Reaction), intent(in) :: this
integer(C_INT) :: sz
sz = reaction_products_size(this % ptr)
end function
function sample_elastic_mu(this, E) result(mu)
class(Reaction), intent(in) :: this
real(C_DOUBLE), intent(in) :: E
real(C_DOUBLE) :: mu
mu = reaction_sample_elastic_mu(this % ptr, E)
end function
function xs(this, temperature, energy) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: temperature
integer(C_INT), intent(in) :: energy
real(C_DOUBLE) :: val
val = reaction_xs(this % ptr, temperature, energy)
end function
function xs_size(this, temperature) result(sz)
class(Reaction), intent(in) :: this
integer(C_INT) :: temperature
integer(C_INT) :: sz
sz = reaction_xs_size(this % ptr, temperature)
end function
function xs_threshold(this, temperature) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: temperature
integer(C_INT) :: val
val = reaction_xs_threshold(this % ptr, temperature)
end function
end module reaction_header

4
src/simulation.F90
View file

@ -2,7 +2,7 @@ module simulation
use, intrinsic :: ISO_C_BINDING
use nuclide_header, only: micro_xs, n_nuclides
use nuclide_header, only: micro_xs, nuclides_size
use photon_header, only: micro_photon_xs, n_elements
implicit none
@ -18,7 +18,7 @@ contains
!$omp parallel
! Allocate array for microscopic cross section cache
allocate(micro_xs(n_nuclides))
allocate(micro_xs(nuclides_size()))
allocate(micro_photon_xs(n_elements))
!$omp end parallel

11
src/simulation.cpp
View file

@ -29,11 +29,7 @@ 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 +72,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 +206,7 @@ int openmc_next_batch(int* status)
initialize_history(&p, simulation::current_work);
// transport particle
transport(&p);
p.transport();
}
// Accumulate time for transport
@ -319,7 +312,7 @@ void initialize_batch()
}
// Reset total starting particle weight used for normalizing tallies
total_weight = 0.0;
simulation::total_weight = 0.0;
// Determine if this batch is the first inactive or active batch.
bool first_inactive = false;

2
src/source.cpp
View file

@ -172,7 +172,7 @@ Bank SourceDistribution::sample() const
if (space_box) {
if (space_box->only_fissionable()) {
// Determine material
auto c = model::cells[cell_index - 1];
auto c = model::cells[cell_index];
auto mat_index = c->material_.size() == 1
? c->material_[0] : c->material_[instance];

81
src/surface.cpp
View file

@ -7,8 +7,9 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/xml_interface.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -340,7 +341,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 +398,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;
@ -1060,12 +1061,11 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
//==============================================================================
extern "C" void
read_surfaces(pugi::xml_node* node)
void read_surfaces(pugi::xml_node node)
{
// Count the number of surfaces.
int n_surfaces = 0;
for (pugi::xml_node surf_node : node->children("surface")) {n_surfaces++;}
for (pugi::xml_node surf_node : node.children("surface")) {n_surfaces++;}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
@ -1075,7 +1075,7 @@ read_surfaces(pugi::xml_node* node)
{
pugi::xml_node surf_node;
int i_surf;
for (surf_node = node->child("surface"), i_surf = 0; surf_node;
for (surf_node = node.child("surface"), i_surf = 0; surf_node;
surf_node = surf_node.next_sibling("surface"), i_surf++) {
std::string surf_type = get_node_value(surf_node, "type", true, true);
@ -1243,59 +1243,30 @@ read_surfaces(pugi::xml_node* node)
}
}
}
// Check to make sure a boundary condition was applied to at least one
// surface
bool boundary_exists = false;
for (const auto& surf : model::surfaces) {
if (surf->bc_ != BC_TRANSMIT) {
boundary_exists = true;
break;
}
}
if (settings::run_mode != RUN_MODE_PLOTTING && !boundary_exists) {
fatal_error("No boundary conditions were applied to any surfaces!");
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind) {return model::surfaces[surf_ind];}
int surface_id(Surface* surf) {return surf->id_;}
int surface_bc(Surface* surf) {return surf->bc_;}
void surface_reflect(Surface* surf, double xyz[3], double uvw[3])
{
Position r {xyz};
Direction u {uvw};
u = surf->reflect(r, u);
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic_;}
bool
surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3],
double uvw[3])
{
Position r {xyz};
Direction u {uvw};
bool rotational = surf->periodic_translate(other, r, u);
// Copy back to arrays
xyz[0] = r.x;
xyz[1] = r.y;
xyz[2] = r.z;
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
return rotational;
}
void free_memory_surfaces_c()
{
for (Surface* surf : model::surfaces) {delete surf;}
model::surfaces.clear();
model::surface_map.clear();
}
int surfaces_size() { return model::surfaces.size(); }
extern "C" void free_memory_surfaces()
{
for (Surface* surf : model::surfaces) {delete surf;}
model::surfaces.clear();
model::surface_map.clear();
}
} // namespace openmc

135
src/surface_header.F90
View file

@ -1,135 +0,0 @@
module surface_header
use, intrinsic :: ISO_C_BINDING
use hdf5_interface
implicit none
interface
pure function surface_pointer(surf_ind) bind(C) result(ptr)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind
type(C_PTR) :: ptr
end function surface_pointer
pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: id
end function surface_id_c
pure function surface_bc_c(surf_ptr) bind(C, name='surface_bc') result(bc)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: bc
end function surface_bc_c
pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) &
bind(C, name='surface_reflect')
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
end subroutine surface_reflect_c
pure function surface_i_periodic_c(surf_ptr) &
bind(C, name="surface_i_periodic") result(i_periodic)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: i_periodic
end function surface_i_periodic_c
function surface_periodic_c(surf_ptr, other_ptr, xyz, uvw) &
bind(C, name="surface_periodic") result(rotational)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
type(C_PTR), intent(in), value :: other_ptr
real(C_DOUBLE), intent(inout) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
logical(C_BOOL) :: rotational
end function surface_periodic_c
subroutine free_memory_surfaces_c() bind(C)
end subroutine free_memory_surfaces_c
end interface
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)
!===============================================================================
type :: Surface
type(C_PTR) :: ptr
contains
procedure :: id => surface_id
procedure :: bc => surface_bc
procedure :: reflect => surface_reflect
procedure :: i_periodic => surface_i_periodic
procedure :: periodic_translate => surface_periodic
end type Surface
type(Surface), allocatable, target :: surfaces(:)
interface
function surfaces_size() result(sz) bind(C)
import C_INT
integer(C_INT) :: sz
end function
end interface
contains
pure function surface_id(this) result(id)
class(Surface), intent(in) :: this
integer(C_INT) :: id
id = surface_id_c(this % ptr)
end function surface_id
pure function surface_bc(this) result(bc)
class(Surface), intent(in) :: this
integer(C_INT) :: bc
bc = surface_bc_c(this % ptr)
end function surface_bc
pure subroutine surface_reflect(this, xyz, uvw)
class(Surface), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
call surface_reflect_c(this % ptr, xyz, uvw)
end subroutine surface_reflect
pure function surface_i_periodic(this) result(i_periodic)
class(Surface), intent(in) :: this
integer(C_INT) :: i_periodic
i_periodic = surface_i_periodic_c(this % ptr)
end function surface_i_periodic
function surface_periodic(this, other, xyz, uvw) result(rotational)
class(Surface), intent(in) :: this
class(Surface), intent(in) :: other
real(C_DOUBLE), intent(inout) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
logical(C_BOOL) :: rotational
rotational = surface_periodic_c(this % ptr, other % ptr, xyz, uvw)
end function surface_periodic
!===============================================================================
! FREE_MEMORY_SURFACES deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_surfaces()
if (allocated(surfaces)) deallocate(surfaces)
call free_memory_surfaces_c()
end subroutine free_memory_surfaces
end module surface_header

18
src/tallies/derivative.cpp
View file

@ -569,9 +569,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.
@ -620,15 +618,7 @@ score_track_derivative(const Particle* p, double distance)
}
}
//! Adjust diff tally flux derivatives for a particle scattering event.
//
//! Note that this subroutine will be called after absorption events in
//! addition to scattering events, but any flux derivatives scored after an
//! 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 +688,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.;
}

75
src/tallies/tally.F90
View file

@ -4,94 +4,19 @@ module tally
use bank_header
use constants
use dict_header, only: EMPTY
use error, only: fatal_error
use geometry_header
use material_header
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

4
src/tallies/tally.cpp
View file

@ -585,9 +585,9 @@ void reduce_tally_results()
// We also need to determine the total starting weight of particles from the
// last realization
double weight_reduced;
MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM,
MPI_Reduce(&simulation::total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM,
0, mpi::intracomm);
if (mpi::master) total_weight = weight_reduced;
if (mpi::master) simulation::total_weight = weight_reduced;
}
#endif

1
src/tallies/tally_header.F90
View file

@ -7,7 +7,6 @@ module tally_header
use dict_header, only: DictIntInt
use hdf5_interface, only: HID_T, HSIZE_T
use message_passing, only: n_procs, master
use nuclide_header, only: nuclide_map_get
use settings, only: reduce_tallies, run_mode
use stl_vector, only: VectorInt
use string, only: to_lower, to_f_string, str_to_int, to_str, to_c_string

55
src/tallies/tally_scoring.cpp
View file

@ -1972,12 +1972,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux,
}
}
//! Score tallies based on a simple count of events (for continuous energy).
//
//! 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]};
@ -2036,12 +2031,7 @@ score_analog_tally_ce(const Particle* p)
match.bins_present_ = false;
}
//! Score tallies based on a simple count of events (for multigroup).
//
//! 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]};
@ -2090,13 +2080,7 @@ score_analog_tally_mg(const Particle* p)
match.bins_present_ = false;
}
//! Score tallies using a tracklength estimate of the flux.
//
//! This is triggered at every event (surface crossing, lattice crossing, or
//! 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
@ -2164,16 +2148,7 @@ score_tracklength_tally(const Particle* p, double distance)
match.bins_present_ = false;
}
//! Score tallies using a 1 / Sigma_t estimate of the flux.
//
//! This is triggered after every collision. It is invalid for tallies that
//! require post-collison information because it can score reactions that didn't
//! actually occur, and we don't a priori know what the outcome will be for
//! 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;
@ -2241,10 +2216,8 @@ score_collision_tally(const Particle* p)
match.bins_present_ = false;
}
//! 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 +2262,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

10
src/track_output.F90
View file

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

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

2
tests/unit_tests/test_capi.py
View file

@ -72,7 +72,7 @@ def test_cell(capi_init):
cell = openmc.capi.cells[1]
assert isinstance(cell.fill, openmc.capi.Material)
cell.fill = openmc.capi.materials[1]
assert str(cell) == 'Cell[1]'
assert str(cell) == 'Cell[0]'
def test_new_cell(capi_init):