crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Convert remainder of tallies to C++

Browse source
This commit is contained in:
Paul Romano 2019-02-21 08:32:32 -06:00
parent 487b62f760
commit 0439eb0fa1
54 changed files with 410 additions and 1420 deletions
Download patch file Download diff file Expand all files Collapse all files

8
CMakeLists.txt
View file

@ -303,7 +303,6 @@ target_compile_options(faddeeva PRIVATE ${cxxflags})
add_library(libopenmc SHARED
src/bank_header.F90
src/api.F90
src/constants.F90
src/dict_header.F90
src/error.F90
@ -321,13 +320,6 @@ add_library(libopenmc SHARED
src/string.F90
src/vector_header.F90
src/xml_interface.F90
src/tallies/tally.F90
src/tallies/tally_derivative_header.F90
src/tallies/tally_filter_header.F90
src/tallies/tally_filter_distribcell.F90
src/tallies/tally_filter_particle.F90
src/tallies/tally_filter_sph_harm.F90
src/tallies/tally_header.F90
src/bank.cpp
src/bremsstrahlung.cpp
src/dagmc.cpp

6
include/openmc/bank.h
View file

@ -42,6 +42,12 @@ extern std::vector<Bank> master_fission_bank;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void free_memory_bank();
} // namespace openmc
#endif // OPENMC_BANK_H

4
include/openmc/capi.h
View file

@ -43,7 +43,7 @@ extern "C" {
int openmc_filter_get_type(int32_t index, const char** type);
int openmc_filter_set_id(int32_t index, int32_t id);
int openmc_finalize();
int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
int openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
void openmc_get_filter_next_id(int32_t* id);
@ -102,7 +102,7 @@ extern "C" {
int openmc_statepoint_write(const char* filename, bool* write_source);
int openmc_tally_allocate(int32_t index, const char* type);
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
int openmc_tally_get_estimator(int32_t index, int* estimator);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n);
int openmc_tally_get_n_realizations(int32_t index, int32_t* n);

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

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

2
include/openmc/cross_sections.h
View file

@ -70,6 +70,8 @@ void read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
//! Read cross_sections.xml and populate data libraries
void read_ce_cross_sections_xml();
void library_clear();
} // namespace openmc
#endif // OPENMC_CROSS_SECTIONS_H

2
include/openmc/dagmc.h
View file

@ -23,7 +23,7 @@ namespace model {
//==============================================================================
void load_dagmc_geometry();
extern "C" void free_memory_dagmc();
void free_memory_dagmc();
void read_geometry_dagmc();
bool read_uwuw_materials(pugi::xml_document& doc);
bool get_uwuw_materials_xml(std::string& s);

2
include/openmc/geometry_aux.h
View file

@ -113,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();
void free_memory_geometry();
} // namespace openmc
#endif // OPENMC_GEOMETRY_AUX_H

2
include/openmc/material.h
View file

@ -112,6 +112,8 @@ private:
//! Read material data from materials.xml
void read_materials_xml();
void free_memory_material();
//==============================================================================
// Fortran compatibility
//==============================================================================

2
include/openmc/mesh.h
View file

@ -130,6 +130,8 @@ void read_meshes(pugi::xml_node root);
//! \param[in] group HDF5 group
void meshes_to_hdf5(hid_t group);
void free_memory_mesh();
} // namespace openmc
#endif // OPENMC_MESH_H

6
include/openmc/nuclide.h
View file

@ -203,6 +203,12 @@ extern MaterialMacroXS material_xs;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void nuclides_clear();
//==============================================================================
// Fortran compatibility
//==============================================================================

2
include/openmc/photon.h
View file

@ -119,6 +119,8 @@ struct ElementMicroXS {
std::pair<double, double> klein_nishina(double alpha);
void free_memory_photon();
//==============================================================================
// Global variables
//==============================================================================

6
include/openmc/settings.h
View file

@ -94,11 +94,15 @@ extern double weight_cutoff; //!< Weight cutoff for Russian roulette
extern double weight_survive; //!< Survival weight after Russian roulette
} // namespace settings
//==============================================================================
// Functions
//==============================================================================
//! Read settings from XML file
//! \param[in] root XML node for <settings>
void read_settings_xml();
extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr);
void free_memory_settings();
} // namespace openmc

5
include/openmc/simulation.h
View file

@ -84,13 +84,12 @@ 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
void free_memory_simulation();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

2
include/openmc/source.h
View file

@ -67,6 +67,8 @@ Bank sample_external_source();
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();
void free_memory_source();
} // namespace openmc
#endif // OPENMC_SOURCE_H

16
include/openmc/surface.h
View file

@ -401,21 +401,7 @@ public:
void read_surfaces(pugi::xml_node node);
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Surface* surface_pointer(int surf_ind);
int surface_id(Surface* surf);
int surface_bc(Surface* surf);
bool surface_sense(Surface* surf, double xyz[3], double uvw[3]);
void surface_reflect(Surface* surf, double xyz[3], double uvw[3]);
int surface_i_periodic(PeriodicSurface* surf);
bool surface_periodic(PeriodicSurface* surf, PeriodicSurface* other,
double xyz[3], double uvw[3]);
void free_memory_surfaces_c();
}
void free_memory_surfaces();
} // namespace openmc
#endif // OPENMC_SURFACE_H

4
include/openmc/tallies/filter.h
View file

@ -102,9 +102,7 @@ namespace model {
Filter* allocate_filter(const std::string& type);
// Filter-related Fortran functions that will be called from C++
extern "C" int verify_filter(int32_t index);
extern "C" Filter* filter_from_f(int32_t index);
extern "C" void filter_update_n_bins(int32_t index);
int verify_filter(int32_t index);
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_H

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

@ -51,6 +51,8 @@ public:
void init_results();
void reset();
void accumulate();
//----------------------------------------------------------------------------
@ -172,7 +174,7 @@ using adaptor_type = xt::xtensor_adaptor<xt::xbuffer_adaptor<double*&, xt::no_ow
void reduce_tally_results();
#endif
extern "C" void free_memory_tally_c();
void free_memory_tally();
} // namespace openmc

2
include/openmc/thermal.h
View file

@ -138,6 +138,8 @@ public:
std::vector<ThermalData> data_;
};
void free_memory_thermal();
} // namespace openmc
#endif // OPENMC_THERMAL_H

14
include/openmc/volume_calc.h
View file

@ -1,5 +1,5 @@
#ifndef VOLUME_CALC_H
#define VOLUME_CALC_H
#ifndef OPENMC_VOLUME_CALC_H
#define OPENMC_VOLUME_CALC_H
#include "openmc/position.h"
@ -66,9 +66,15 @@ private:
//==============================================================================
namespace model {
extern std::vector<VolumeCalculation> volume_calcs;
extern std::vector<VolumeCalculation> volume_calcs;
}
//==============================================================================
// Non-member functions
//==============================================================================
void free_memory_volume();
} // namespace openmc
#endif // VOLUME_CALC_H
#endif // OPENMC_VOLUME_CALC_H

7
openmc/capi/tally.py
View file

@ -32,7 +32,7 @@ _dll.openmc_tally_allocate.errcheck = _error_handler
_dll.openmc_tally_get_active.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_active.restype = c_int
_dll.openmc_tally_get_active.errcheck = _error_handler
_dll.openmc_tally_get_estimator.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_estimator.argtypes = [c_int32, POINTER(int32)]
_dll.openmc_tally_get_estimator.restype = c_int
_dll.openmc_tally_get_estimator.errcheck = _error_handler
_dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
@ -84,6 +84,7 @@ _dll.openmc_tally_set_scores.errcheck = _error_handler
_dll.openmc_tally_set_type.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_set_type.restype = c_int
_dll.openmc_tally_set_type.errcheck = _error_handler
_dll.tallies_size.restype = c_size_t
_SCORES = {
@ -227,7 +228,7 @@ class Tally(_FortranObjectWithID):
@property
def estimator(self):
estimator = c_int32()
estimator = c_int()
_dll.openmc_tally_get_estimator(self._index, estimator)
return _ESTIMATORS[estimator.value]
@ -387,7 +388,7 @@ class _TallyMapping(Mapping):
yield Tally(index=i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_tallies').value
return _dll.tallies_size()
def __repr__(self):
return repr(dict(self))

136
src/api.F90
View file

@ -1,136 +0,0 @@
module openmc_api
use, intrinsic :: ISO_C_BINDING
use constants
use error
use geometry, only: find_cell
use particle_header
use string, only: to_str
#ifdef DAGMC
use dagmc_header, only: free_memory_dagmc
#endif
implicit none
private
contains
!===============================================================================
! OPENMC_FIND_CELL determines what cell contains a given point in space
!===============================================================================
function openmc_find_cell(xyz, index, instance) result(err) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3) ! Cartesian point
integer(C_INT32_T), intent(out) :: index
integer(C_INT32_T), intent(out) :: instance
integer(C_INT) :: err
logical :: found
type(Particle) :: p
call particle_initialize(p)
p % coord(1) % xyz(:) = xyz
p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
call find_cell(p, found)
index = -1
instance = -1
err = E_UNASSIGNED
if (found) then
index = p % coord(p % n_coord) % cell
instance = p % cell_instance
err = 0
else
err = E_GEOMETRY
call set_errmsg("Could not find cell at position (" // &
trim(to_str(xyz(1))) // "," // trim(to_str(xyz(2))) // "," // &
trim(to_str(xyz(3))) // ").")
end if
end function openmc_find_cell
!===============================================================================
! FREE_MEMORY deallocates and clears all global allocatable arrays in the
! program
!===============================================================================
subroutine free_memory() bind(C)
use bank_header
use material_header
use settings
use simulation_header
use tally_filter_header
use tally_header
interface
subroutine free_memory_source() bind(C)
end subroutine
subroutine free_memory_material() bind(C)
end subroutine
subroutine free_memory_mesh() bind(C)
end subroutine free_memory_mesh
subroutine free_memory_settings() bind(C)
end subroutine free_memory_settings
subroutine free_memory_bank() bind(C)
end subroutine free_memory_bank
subroutine free_memory_cmfd() bind(C)
end subroutine free_memory_cmfd
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 free_memory_photon() bind(C)
end subroutine
subroutine sab_clear() bind(C)
end subroutine
subroutine library_clear() bind(C)
end subroutine
subroutine nuclides_clear() bind(C)
end subroutine
#ifdef DAGMC
subroutine free_memory_dagmc() bind(C)
end subroutine
#endif
end interface
call free_memory_geometry()
call free_memory_surfaces()
call free_memory_material()
call free_memory_volume()
call free_memory_simulation()
call free_memory_photon()
call free_memory_settings()
call sab_clear()
call library_clear()
call nuclides_clear()
call free_memory_source()
call free_memory_mesh()
call free_memory_tally()
call free_memory_tally_filter()
call free_memory_bank()
call free_memory_cmfd()
#ifdef DAGMC
call free_memory_dagmc()
#endif
end subroutine free_memory
end module openmc_api

28
src/bank.cpp
View file

@ -26,6 +26,22 @@ std::vector<Bank> master_fission_bank;
} // namespace simulation
//==============================================================================
// Non-member functions
//==============================================================================
void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
{
simulation::fission_bank.clear();
}
#ifdef _OPENMP
simulation::master_fission_bank.clear();
#endif
}
//==============================================================================
// C API
//==============================================================================
@ -58,18 +74,6 @@ extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n)
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_bank()
{
simulation::source_bank.clear();
#pragma omp parallel
{
simulation::fission_bank.clear();
}
#ifdef _OPENMP
simulation::master_fission_bank.clear();
#endif
}
extern "C" int fission_bank_delayed_group(int64_t i) {
return simulation::fission_bank[i-1].delayed_group;
}

11
src/cmfd_solver.cpp
View file

@ -1,3 +1,5 @@
#include "openmc/cmfd_solver.h"
#include <vector>
#include <cmath>
@ -348,11 +350,7 @@ int openmc_run_linsolver(const double* A_data, const double* b, double* x,
}
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_cmfd()
void free_memory_cmfd()
{
cmfd::indptr.clear();
cmfd::indices.clear();
@ -360,5 +358,4 @@ extern "C" void free_memory_cmfd()
cmfd::indexmap.resize({0});
}
} // namespace openmc
} // namespace openmc

6
src/cross_sections.cpp
View file

@ -417,11 +417,7 @@ void read_ce_cross_sections_xml()
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void library_clear() {
void library_clear() {
data::libraries.clear();
data::library_map.clear();
}

46
src/finalize.cpp
View file

@ -1,23 +1,58 @@
#include "openmc/finalize.h"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/cmfd_solver.h"
#include "openmc/constants.h"
#include "openmc/cross_sections.h"
#include "openmc/dagmc.h"
#include "openmc/eigenvalue.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/material.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/surface.h"
#include "openmc/thermal.h"
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
#include "openmc/volume_calc.h"
#include "xtensor/xview.hpp"
using namespace openmc;
namespace openmc {
// Functions defined in Fortran
extern "C" void free_memory();
void free_memory()
{
free_memory_geometry();
free_memory_surfaces();
free_memory_material();
free_memory_volume();
free_memory_simulation();
free_memory_photon();
free_memory_settings();
free_memory_thermal();
library_clear();
nuclides_clear();
free_memory_source();
free_memory_mesh();
free_memory_tally();
free_memory_bank();
free_memory_cmfd();
#ifdef DAGMC
free_memory_dagmc();
#endif
}
}
using namespace openmc;
int openmc_finalize()
{
@ -81,7 +116,6 @@ int openmc_finalize()
data::energy_max = {INFTY, INFTY};
data::energy_min = {0.0, 0.0};
n_tallies = 0;
model::root_universe = -1;
openmc_set_seed(DEFAULT_SEED);
@ -98,8 +132,8 @@ int openmc_finalize()
int openmc_reset()
{
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_reset(i);
for (auto& t : model::tallies) {
t->reset();
}
// Reset global tallies

28
src/geometry.cpp
View file

@ -464,4 +464,32 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed,
}
}
//==============================================================================
// C API
//==============================================================================
extern "C" int
openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance)
{
Particle p;
p.initialize();
std::copy(xyz, xyz + 3, p.coord[0].xyz);
p.coord[0].uvw[0] = 0.0;
p.coord[0].uvw[1] = 0.0;
p.coord[0].uvw[2] = 1.0;
if (!find_cell(&p, false)) {
std::stringstream msg;
msg << "Could not find cell at position (" << xyz[0] << ", " << xyz[1]
<< ", " << xyz[2] << ").";
set_errmsg(msg);
return OPENMC_E_GEOMETRY;
}
*index = p.coord[p.n_coord-1].cell;
*instance = p.cell_instance;
return 0;
}
} // namespace openmc

15
src/material.cpp
View file

@ -913,6 +913,13 @@ void read_materials_xml()
}
}
void free_memory_material()
{
for (Material *mat : model::materials) {delete mat;}
model::materials.clear();
model::material_map.clear();
}
//==============================================================================
// C API
//==============================================================================
@ -1123,6 +1130,7 @@ openmc_material_set_volume(int32_t index, double volume)
extern "C" int
openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end)
{
// TODO: off-by-one
if (index_start) *index_start = model::materials.size() + 1;
if (index_end) *index_end = model::materials.size() + n;
for (int32_t i = 0; i < n; i++) {
@ -1159,13 +1167,6 @@ extern "C" {
{
return model::materials[i_mat - 1]->density_gpcc_;
}
void free_memory_material()
{
for (Material *mat : model::materials) {delete mat;}
model::materials.clear();
model::material_map.clear();
}
}
} // namespace openmc

11
src/mesh.cpp
View file

@ -894,6 +894,12 @@ void meshes_to_hdf5(hid_t group)
close_group(meshes_group);
}
void free_memory_mesh()
{
model::meshes.clear();
model::mesh_map.clear();
}
//==============================================================================
// Fortran compatibility
//==============================================================================
@ -937,11 +943,6 @@ extern "C" {
m->get_indices_from_bin(bin, ijk);
}
void free_memory_mesh()
{
model::meshes.clear();
model::mesh_map.clear();
}
}

12
src/nuclide.cpp
View file

@ -956,6 +956,12 @@ openmc_nuclide_name(int index, const char** name)
}
}
void nuclides_clear()
{
data::nuclides.clear();
data::nuclide_map.clear();
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -966,10 +972,4 @@ extern "C" bool multipole_in_range(const Nuclide* nuc, double E)
E <= nuc->multipole_->E_max_;
}
extern "C" void nuclides_clear()
{
data::nuclides.clear();
data::nuclide_map.clear();
}
} // namespace openmc

6
src/photon.cpp
View file

@ -771,11 +771,7 @@ std::pair<double, double> klein_nishina(double alpha)
return {alpha_out, mu};
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_photon()
void free_memory_photon()
{
data::elements.clear();
data::compton_profile_pz.resize({0});

11
src/settings.cpp
View file

@ -782,6 +782,12 @@ void read_settings_xml()
}
}
void free_memory_settings() {
settings::statepoint_batch.clear();
settings::sourcepoint_batch.clear();
settings::res_scat_nuclides.clear();
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -794,11 +800,6 @@ extern "C" {
return settings::path_input.c_str();
}
void free_memory_settings() {
settings::statepoint_batch.clear();
settings::sourcepoint_batch.clear();
settings::res_scat_nuclides.clear();
}
}
} // namespace openmc

19
src/simulation.cpp
View file

@ -163,8 +163,8 @@ int openmc_simulation_finalize()
}
// Deactivate all tallies
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_set_active(i, false);
for (auto& t : model::tallies) {
t->active_ = false;
}
// Stop timers and show timing statistics
@ -342,9 +342,8 @@ void initialize_batch()
} else if (first_active) {
simulation::time_inactive.stop();
simulation::time_active.start();
for (int i = 1; i <= n_tallies; ++i) {
// TODO: change one-based index
openmc_tally_set_active(i, true);
for (auto& t : model::tallies) {
t->active_ = true;
}
}
@ -544,11 +543,11 @@ void calculate_work()
#ifdef OPENMC_MPI
void broadcast_results() {
// Broadcast tally results so that each process has access to results
for (int i = 0; i < n_tallies; ++i) {
for (auto& t : model::tallies) {
// Create a new datatype that consists of all values for a given filter
// bin and then use that to broadcast. This is done to minimize the
// chance of the 'count' argument of MPI_BCAST exceeding 2**31
auto& results = model::tallies[i]->results_;
auto& results = t->results_;
auto shape = results.shape();
int count_per_filter = shape[1] * shape[2];
@ -575,6 +574,12 @@ void broadcast_results() {
#endif
void free_memory_simulation()
{
simulation::k_generation.clear();
simulation::entropy.clear();
}
//==============================================================================
// Fortran compatibility
//==============================================================================

10
src/source.cpp
View file

@ -327,15 +327,15 @@ Bank sample_external_source()
return site;
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void free_memory_source()
void free_memory_source()
{
model::external_sources.clear();
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {

9
src/state_point.cpp
View file

@ -235,9 +235,8 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_attribute(file_id, "tallies_present", 1);
// Write all tally results
for (int i = 0; i < model::tallies.size(); ++i) {
for (const auto& tally : model::tallies) {
// Write sum and sum_sq for each bin
const auto& tally {model::tallies[i]};
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
auto& results = tally->results_;
@ -426,9 +425,8 @@ void load_state_point()
if (present) {
hid_t tallies_group = open_group(file_id, "tallies");
for (int i = 0; i < model::tallies.size(); ++i) {
for (auto& tally : model::tallies) {
// Read sum, sum_sq, and N for each bin
auto& tally {model::tallies[i]};
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
auto& results = tally->results_;
@ -697,9 +695,8 @@ void write_tally_results_nr(hid_t file_id)
write_dataset(file_id, "global_tallies", gt);
}
for (int i = 0; i < n_tallies; ++i) {
for (const auto& t : model::tallies) {
// Skip any tallies that are not active
const auto& t {model::tallies[i]};
if (!t->active_) continue;
if (mpi::master && !object_exists(file_id, "tallies_present")) {

6
src/surface.cpp
View file

@ -1258,11 +1258,7 @@ void read_surfaces(pugi::xml_node node)
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void free_memory_surfaces()
void free_memory_surfaces()
{
for (Surface* surf : model::surfaces) {delete surf;}
model::surfaces.clear();

62
src/tallies/filter.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/tallies/filter.h"
#include <algorithm> // for max
#include <string>
#include "openmc/capi.h"
@ -105,6 +106,43 @@ allocate_filter(const std::string& type)
// C API functions
//==============================================================================
int verify_filter(int32_t index)
{
if (index < 1 || index > model::tally_filters.size()) {
set_errmsg("Filter index is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
extern "C" int
openmc_filter_get_id(int32_t index, int32_t* id)
{
int err = verify_filter(index);
if (err) return err;
// TODO: off-by-one
*id = model::tally_filters[index-1]->id_;
return 0;
}
extern "C" int
openmc_filter_set_id(int32_t index, int32_t id)
{
int err = verify_filter(index);
if (err) return err;
if (model::filter_map.find(id) != model::filter_map.end()) {
set_errmsg("Two filters have the same ID: " + std::to_string(id));
return OPENMC_E_INVALID_ID;
}
// TODO: off-by-one
model::tally_filters[index-1]->id_ = id;
model::filter_map[id] = index - 1;
return 0;
}
extern "C" int
openmc_filter_get_type(int32_t index, const char** type)
{
@ -115,6 +153,30 @@ openmc_filter_get_type(int32_t index, const char** type)
*type = model::tally_filters[index-1]->type().c_str();
}
extern "C" int
openmc_get_filter_index(int32_t id, int32_t* index)
{
auto it = model::filter_map.find(id);
if (it == model::filter_map.end()) {
set_errmsg("No filter exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
}
// TODO: off-by-one
*index = it->second + 1;
return 0;
}
extern "C" void
openmc_get_filter_next_id(int32_t* id)
{
int32_t largest_filter_id = 0;
for (const auto& t : model::tally_filters) {
largest_filter_id = std::max(largest_filter_id, t->id_);
}
*id = largest_filter_id + 1;
}
extern "C" int
openmc_new_filter(const char* type, int32_t* index)
{

3
src/tallies/filter_cell.cpp
View file

@ -77,7 +77,8 @@ openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n)
int err = verify_filter(index);
if (err) return err;
auto filt = filter_from_f(index);
// TODO: off-by-one
const auto& filt = model::tally_filters[index-1].get();
if (filt->type() != "cell") {
set_errmsg("Tried to get cells from a non-cell filter.");
return OPENMC_E_INVALID_TYPE;

5
src/tallies/filter_energy.cpp
View file

@ -124,7 +124,7 @@ openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<EnergyFilter*>(filt_base);
// Check the filter type.
@ -147,7 +147,7 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<EnergyFilter*>(filt_base);
// Check the filter type.
@ -161,7 +161,6 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
filt->bins_.resize(n);
for (int i = 0; i < n; i++) filt->bins_[i] = energies[i];
filt->n_bins_ = n - 1;
filter_update_n_bins(index);
return 0;
}

5
src/tallies/filter_legendre.cpp
View file

@ -53,7 +53,7 @@ openmc_legendre_filter_get_order(int32_t index, int* order)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<LegendreFilter*>(filt_base);
// Check the filter type.
@ -75,7 +75,7 @@ openmc_legendre_filter_set_order(int32_t index, int order)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<LegendreFilter*>(filt_base);
// Check the filter type.
@ -87,7 +87,6 @@ openmc_legendre_filter_set_order(int32_t index, int order)
// Update the filter.
filt->order_ = order;
filt->n_bins_ = order + 1;
filter_update_n_bins(index);
return 0;
}

5
src/tallies/filter_material.cpp
View file

@ -78,7 +78,7 @@ openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<MaterialFilter*>(filt_base);
// Check the filter type.
@ -101,7 +101,7 @@ openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<MaterialFilter*>(filt_base);
// Check the filter type.
@ -117,7 +117,6 @@ openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins)
filt->n_bins_ = filt->materials_.size();
filt->map_.clear();
for (int i = 0; i < n; i++) filt->map_[filt->materials_[i]] = i;
filter_update_n_bins(index);
return 0;
}

5
src/tallies/filter_mesh.cpp
View file

@ -91,7 +91,7 @@ openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<MeshFilter*>(filt_base);
// Check the filter type.
@ -113,7 +113,7 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
if (err) return err;
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<MeshFilter*>(filt_base);
// Check the filter type.
@ -130,7 +130,6 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
// Update the filter.
filt->set_mesh(index_mesh);
filter_update_n_bins(index);
return 0;
}

3
src/tallies/filter_sph_harm.cpp
View file

@ -100,7 +100,7 @@ check_sphharm_filter(int32_t index)
}
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<SphericalHarmonicsFilter*>(filt_base);
// Check the filter type.
@ -155,7 +155,6 @@ openmc_sphharm_filter_set_order(int32_t index, int order)
// Update the filter.
filt->order_ = order;
filt->n_bins_ = (order + 1) * (order + 1);
filter_update_n_bins(index);
return 0;
}

3
src/tallies/filter_sptl_legendre.cpp
View file

@ -107,7 +107,7 @@ check_sptl_legendre_filter(int32_t index)
}
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<SpatialLegendreFilter*>(filt_base);
// Check the filter type.
@ -161,7 +161,6 @@ openmc_spatial_legendre_filter_set_order(int32_t index, int order)
// Update the filter.
filt->order_ = order;
filt->n_bins_ = order + 1;
filter_update_n_bins(index);
return 0;
}

3
src/tallies/filter_zernike.cpp
View file

@ -132,7 +132,7 @@ check_zernike_filter(int32_t index)
}
// Get a pointer to the filter and downcast.
auto* filt_base = filter_from_f(index);
const auto& filt_base = model::tally_filters[index-1].get();
auto* filt = dynamic_cast<ZernikeFilter*>(filt_base);
// Check the filter type.
@ -185,7 +185,6 @@ openmc_zernike_filter_set_order(int32_t index, int order)
// Update the filter.
filt->set_order(order);
filter_update_n_bins(index);
return 0;
}

69
src/tallies/tally.F90
View file

@ -1,69 +0,0 @@
module tally
use, intrinsic :: ISO_C_BINDING
use bank_header
use constants
use material_header
use message_passing
use settings
use simulation_header
use tally_derivative_header
use tally_filter
use tally_header
implicit none
contains
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_tally_allocate(index, type) result(err) bind(C)
! Set the type of the tally
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR), intent(in) :: type(*)
integer(C_INT) :: err
integer(C_INT32_T) :: empty(0)
character(:), allocatable :: type_
interface
function tally_pointer(indx) bind(C) result(ptr)
import C_INT, C_PTR
integer(C_INT), value :: indx
type(C_PTR) :: ptr
end function
end interface
! Convert C string to Fortran string
type_ = to_f_string(type)
err = 0
if (index >= 1 .and. index <= n_tallies) then
if (allocated(tallies(index) % obj)) then
err = E_ALLOCATE
call set_errmsg("Tally type has already been set.")
else
select case (type_)
case ('generic')
allocate(TallyObject :: tallies(index) % obj)
case default
err = E_UNASSIGNED
call set_errmsg("Unknown tally type: " // trim(type_))
end select
! Assign the pointer to the C++ tally
tallies(index) % obj % ptr = tally_pointer(index - 1)
! When a tally is allocated, set it to have 0 filters
err = openmc_tally_set_filters(index, 0, empty)
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in tallies array is out of bounds.")
end if
end function openmc_tally_allocate
end module tally

138
src/tallies/tally.cpp
View file

@ -31,8 +31,9 @@
#include "xtensor/xbuilder.hpp" // for empty_like
#include "xtensor/xview.hpp"
#include <algorithm> // for max
#include <array>
#include <cstddef>
#include <cstddef> // for size_t
#include <sstream>
#include <string>
@ -519,6 +520,15 @@ void Tally::init_results()
results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
}
void Tally::reset()
{
n_realizations_ = 0;
// TODO: Change to zero when xtensor is updated
if (results_.size() != 1) {
xt::view(results_, xt::all()) = 0.0;
}
}
void Tally::accumulate()
{
// Increment number of realizations
@ -1004,8 +1014,8 @@ setup_active_tallies()
}
}
extern "C" void
free_memory_tally_c()
void
free_memory_tally()
{
#pragma omp parallel
{
@ -1013,6 +1023,7 @@ free_memory_tally_c()
}
model::tally_filters.clear();
model::filter_map.clear();
model::tallies.clear();
@ -1022,12 +1033,119 @@ free_memory_tally_c()
model::active_collision_tallies.clear();
model::active_meshsurf_tallies.clear();
model::active_surface_tallies.clear();
model::tally_map.clear();
}
//==============================================================================
// C-API functions
//==============================================================================
extern "C" int
openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end)
{
// TODO: off-by-one
if (index_start) *index_start = model::tallies.size() + 1;
if (index_end) *index_end = model::tallies.size() + n;
for (int i = 0; i < n; ++i) {
model::tallies.push_back(std::make_unique<Tally>());
}
return 0;
}
extern "C" int
openmc_get_tally_index(int32_t id, int32_t* index)
{
auto it = model::tally_map.find(id);
if (it == model::tally_map.end()) {
set_errmsg("No tally exists with ID=" + std::to_string(id) + ".");
return OPENMC_E_INVALID_ID;
}
*index = it->second;
return 0;
}
extern "C" void
openmc_get_tally_next_id(int32_t* id)
{
int32_t largest_tally_id = 0;
for (const auto& t : model::tallies) {
largest_tally_id = std::max(largest_tally_id, t->id_);
}
*id = largest_tally_id + 1;
}
extern "C" int
openmc_tally_get_estimator(int32_t index, int* estimator)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
*estimator = model::tallies[index-1]->estimator_;
return 0;
}
extern "C" int
openmc_tally_set_estimator(int32_t index, const char* estimator)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
auto& t {model::tallies[index-1]};
std::string est = estimator;
if (est == "analog") {
t->estimator_ = ESTIMATOR_ANALOG;
} else if (est == "collision") {
t->estimator_ = ESTIMATOR_COLLISION;
} else if (est == "tracklength") {
t->estimator_ = ESTIMATOR_TRACKLENGTH;
} else {
set_errmsg("Unknown tally estimator: " + est);
return OPENMC_E_INVALID_ARGUMENT;
}
return 0;
}
extern "C" int
openmc_tally_get_id(int32_t index, int32_t* id)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
*id = model::tallies[index-1]->id_;
return 0;
}
extern "C" int
openmc_tally_set_id(int32_t index, int32_t id)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
if (model::tally_map.find(id) != model::tally_map.end()) {
set_errmsg("Two or more tallies use the same unique ID: "
+ std::to_string(id));
return OPENMC_E_INVALID_ID;
}
// TODO: off-by-one
model::tallies[index-1]->id_ = id;
return 0;
}
extern "C" int
openmc_tally_get_type(int32_t index, int32_t* type)
{
@ -1216,12 +1334,7 @@ openmc_tally_reset(int32_t index)
}
// TODO: off-by-one
auto& t {model::tallies[index-1]};
t->n_realizations_ = 0;
// TODO: Change to zero when xtensor is updated
if (t->results_.size() != 1) {
xt::view(t->results_, xt::all()) = 0.0;
}
model::tallies[index-1]->reset();
return 0;
}
@ -1283,13 +1396,6 @@ extern "C" size_t tallies_size() { return model::tallies.size(); }
extern "C" {
Tally* tally_pointer(int indx) {return model::tallies[indx].get();}
void
extend_tallies_c(int n)
{
for (int i = 0; i < n; ++i)
model::tallies.push_back(std::make_unique<Tally>());
}
int active_tallies_data(int i)
{return model::active_tallies[i-1];}

33
src/tallies/tally_derivative_header.F90
View file

@ -1,33 +0,0 @@
module tally_derivative_header
use, intrinsic :: ISO_C_BINDING
implicit none
!===============================================================================
! TALLYDERIVATIVE describes a first-order derivative that can be applied to
! tallies.
!===============================================================================
type, bind(C), public :: TallyDerivative
integer(C_INT) :: id
integer(C_INT) :: variable
integer(C_INT) :: diff_material
integer(C_INT) :: diff_nuclide
real(C_DOUBLE) :: flux_deriv
end type TallyDerivative
interface
function n_tally_derivs() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
function tally_deriv_c(i) result(deriv) bind(C)
import C_INT, TallyDerivative
integer(C_INT), value :: i
type(TallyDerivative), pointer :: deriv
end function
end interface
end module tally_derivative_header

21
src/tallies/tally_filter_distribcell.F90
View file

@ -1,21 +0,0 @@
module tally_filter_distribcell
use tally_filter_header
implicit none
private
type, public, extends(TallyFilter) :: DistribcellFilter
contains
procedure :: initialize => initialize_distribcell
end type DistribcellFilter
contains
subroutine initialize_distribcell(this)
class(DistribcellFilter), intent(inout) :: this
call this % initialize_cpp()
this % n_bins = this % n_bins_cpp()
end subroutine initialize_distribcell
end module tally_filter_distribcell

370
src/tallies/tally_filter_header.F90
View file

@ -1,370 +0,0 @@
module tally_filter_header
use, intrinsic :: ISO_C_BINDING
use constants, only: MAX_LINE_LEN
use dict_header, only: DictIntInt
use error
use hdf5_interface, only: HID_T
use particle_header, only: Particle
use string, only: to_str
use xml_interface, only: XMLNode
implicit none
private
public :: free_memory_tally_filter
public :: verify_filter
public :: openmc_extend_filters
public :: openmc_filter_get_id
public :: openmc_filter_set_id
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
! example, a cell filter indicates that only particles in a specified cell
! should score to the tally.
!===============================================================================
type, public, abstract :: TallyFilter
integer :: n_bins = 0
type(C_PTR) :: ptr
contains
procedure :: id => filter_get_id_f90
procedure :: from_xml
procedure :: to_statepoint
procedure :: initialize
procedure :: n_bins_cpp
procedure :: from_xml_cpp
procedure :: initialize_cpp
end type TallyFilter
!===============================================================================
! Pure C++ filters
!===============================================================================
type, public, extends(TallyFilter) :: AzimuthalFilter
end type
type, public, extends(TallyFilter) :: CellFilter
end type
type, public, extends(TallyFilter) :: CellbornFilter
end type
type, public, extends(CellFilter) :: CellFromFilter
end type
type, public, extends(TallyFilter) :: DelayedGroupFilter
end type
type, public, extends(TallyFilter) :: EnergyFilter
end type
type, public, extends(EnergyFilter) :: EnergyoutFilter
end type
type, public, extends(TallyFilter) :: EnergyFunctionFilter
end type
type, public, extends(TallyFilter) :: LegendreFilter
end type
type, public, extends(TallyFilter) :: MaterialFilter
end type
type, public, extends(TallyFilter) :: MeshFilter
end type
type, public, extends(TallyFilter) :: MeshSurfaceFilter
end type
type, public, extends(TallyFilter) :: MuFilter
end type
type, public, extends(TallyFilter) :: PolarFilter
end type
type, public, extends(TallyFilter) :: SpatialLegendreFilter
end type
type, public, extends(TallyFilter) :: SurfaceFilter
! True if this filter is used for surface currents
logical :: current = .false.
end type
type, public, extends(TallyFilter) :: UniverseFilter
end type
type, public, extends(TallyFilter) :: ZernikeFilter
end type
type, public, extends(ZernikeFilter) :: ZernikeRadialFilter
end type
!===============================================================================
! TALLYFILTERCONTAINER contains an allocatable TallyFilter object for arrays of
! TallyFilters
!===============================================================================
type, public :: TallyFilterContainer
class(TallyFilter), allocatable :: obj
end type TallyFilterContainer
integer(C_INT32_T), public, bind(C) :: n_filters = 0 ! # of filters
type(TallyFilterContainer), public, allocatable, target :: filters(:)
! Dictionary that maps user IDs to indices in 'filters'
type(DictIntInt), public :: filter_dict
! The largest filter ID that has been specified in the system. This is useful
! in case the code needs to find an ID for a new filter.
integer :: largest_filter_id
contains
!===============================================================================
! TallyFilter implementation
!===============================================================================
function filter_get_id_f90(this) result(id)
class(TallyFilter) :: this
integer(C_INT32_T) :: id
interface
function filter_get_id(filt) result(id) bind(C)
import C_PTR, C_INT32_T
type(C_PTR), value :: filt
integer(C_INT32_T) :: id
end function
end interface
id = filter_get_id(this % ptr)
end function
subroutine from_xml(this, node)
class(TallyFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
call this % from_xml_cpp(node)
this % n_bins = this % n_bins_cpp()
end subroutine from_xml
subroutine to_statepoint(this, filter_group)
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
interface
subroutine filter_to_statepoint(filt, filter_group) bind(C)
import C_PTR, HID_T
type(C_PTR), value :: filt
integer(HID_T), intent(in), value :: filter_group
end subroutine filter_to_statepoint
end interface
call filter_to_statepoint(this % ptr, filter_group)
end subroutine to_statepoint
subroutine initialize(this)
class(TallyFilter), intent(inout) :: this
call this % initialize_cpp()
end subroutine initialize
function n_bins_cpp(this) result(n_bins)
class(TallyFilter), intent(in) :: this
integer :: n_bins
interface
function filter_n_bins(filt) result(n_bins) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: n_bins
end function filter_n_bins
end interface
n_bins = filter_n_bins(this % ptr)
end function n_bins_cpp
subroutine from_xml_cpp(this, node)
class(TallyFilter), intent(inout) :: this
class(XMLNode), intent(in) :: node
interface
subroutine filter_from_xml(filt, node) bind(C)
import C_PTR
type(C_PTR), value :: filt
type(C_PTR) :: node
end subroutine filter_from_xml
end interface
call filter_from_xml(this % ptr, node % ptr)
end subroutine from_xml_cpp
subroutine initialize_cpp(this)
class(TallyFilter), intent(inout) :: this
interface
subroutine filter_initialize(filt) bind(C)
import C_PTR
type(C_PTR), value :: filt
end subroutine filter_initialize
end interface
call filter_initialize(this % ptr)
end subroutine initialize_cpp
!===============================================================================
! FILTER_FROM_F given a Fortran index, return a pointer to a C++ filter.
!===============================================================================
function filter_from_f(index) result(filt) bind(C)
integer(C_INT32_T), intent(in), value :: index
type(C_PTR) :: filt
filt = filters(index) % obj % ptr
end function
!===============================================================================
! FILTER_UPDATE_N_BINS given a Fortran index, updates filt % n_bins using C++.
!===============================================================================
subroutine filter_update_n_bins(index) bind(C)
integer(C_INT32_T), intent(in), value :: index
filters(index) % obj % n_bins = filters(index) % obj % n_bins_cpp()
end subroutine
!===============================================================================
! FREE_MEMORY_TALLY_FILTER deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_tally_filter()
n_filters = 0
if (allocated(filters)) deallocate(filters)
call filter_dict % clear()
largest_filter_id = 0
end subroutine free_memory_tally_filter
!===============================================================================
! VERIFY_FILTER makes sure that given a filter index, the size of the filters
! array is sufficient and a filter object has already been allocated.
!===============================================================================
function verify_filter(index) result(err) bind(C)
integer(C_INT32_T), intent(in), value :: index
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= n_filters) then
if (.not. allocated(filters(index) % obj)) then
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
end if
end function verify_filter
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_filters(n, index_start, index_end) result(err) bind(C)
! Creates or extends the filters array
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
integer :: i ! loop counter
type(TallyFilterContainer), allocatable :: temp(:) ! temporary filters
if (n_filters == 0) then
! Allocate filters array
allocate(filters(n))
else
! Move filters to temporary array
allocate(temp(n_filters + n))
do i = 1, n_filters
call move_alloc(filters(i) % obj, temp(i) % obj)
end do
! Move filters back from temporary array to filters array
call move_alloc(temp, filters)
end if
! Return indices in filters array
if (present(index_start)) index_start = n_filters + 1
if (present(index_end)) index_end = n_filters + n
n_filters = n_filters + n
err = 0
end function openmc_extend_filters
function openmc_filter_get_id(index, id) result(err) bind(C)
! Return the ID of a filter
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= n_filters) then
id = filters(index) % obj % id()
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
end if
end function openmc_filter_get_id
function openmc_filter_set_id(index, id) result(err) bind(C)
! Set the ID of a filter
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
interface
subroutine filter_set_id(filt, id) bind(C)
import C_PTR, C_INT32_T
type(C_PTR), value :: filt
integer(C_INT32_T), value :: id
end subroutine
end interface
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
call filter_set_id(filters(index) % obj % ptr, id)
call filter_dict % set(id, index)
if (id > largest_filter_id) largest_filter_id = id
err = 0
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
end if
end function openmc_filter_set_id
function openmc_get_filter_index(id, index) result(err) bind(C)
! Returns the index in the filters array of a filter with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(filters)) then
if (filter_dict % has(id)) then
index = filter_dict % get(id)
err = 0
else
err = E_INVALID_ID
call set_errmsg("No filter exists with ID=" // trim(to_str(id)) // ".")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory has not been allocated for filters.")
end if
end function openmc_get_filter_index
subroutine openmc_get_filter_next_id(id) bind(C)
! Returns an ID number that has not been used by any other filters.
integer(C_INT32_T), intent(out) :: id
id = largest_filter_id + 1
end subroutine openmc_get_filter_next_id
end module tally_filter_header

36
src/tallies/tally_filter_particle.F90
View file

@ -1,36 +0,0 @@
module tally_filter_particle
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
implicit none
private
!===============================================================================
! PARTICLEFILTER specifies what particles can score to a tally
!===============================================================================
type, public, extends(TallyFilter) :: ParticleFilter
contains
procedure :: particles
end type ParticleFilter
contains
function particles(this, i) result(ptype)
class(ParticleFilter), intent(in) :: this
integer, intent(in) :: i
integer :: ptype
interface
function particle_filter_particles(filt, i) result(ptype) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT), value :: i
integer(C_INT) :: ptype
end function
end interface
ptype = particle_filter_particles(this % ptr, i)
end function particles
end module tally_filter_particle

33
src/tallies/tally_filter_sph_harm.F90
View file

@ -1,33 +0,0 @@
module tally_filter_sph_harm
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
implicit none
private
integer, public, parameter :: COSINE_SCATTER = 1
integer, public, parameter :: COSINE_PARTICLE = 2
type, public, extends(TallyFilter) :: SphericalHarmonicsFilter
contains
procedure :: cosine
end type SphericalHarmonicsFilter
contains
function cosine(this) result(val)
class(SphericalHarmonicsFilter) :: this
integer :: val
interface
function sphharm_filter_get_cosine(filt) result(val) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: val
end function
end interface
val = sphharm_filter_get_cosine(this % ptr)
end function cosine
end module tally_filter_sph_harm

551
src/tallies/tally_header.F90
View file

@ -1,551 +0,0 @@
module tally_header
use, intrinsic :: ISO_C_BINDING
use constants
use error
use dict_header, only: DictIntInt
use hdf5_interface, only: HID_T, HSIZE_T
use message_passing, only: n_procs, master
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
use tally_filter_header, only: TallyFilterContainer, filters, n_filters
use tally_filter
implicit none
interface
function openmc_tally_set_filters(index, n, filter_indices) 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_INT32_T), intent(in) :: filter_indices(n)
integer(C_INT) :: err
end function
function openmc_tally_set_active(index, active) result(err) bind(C)
import C_INT32_T, C_BOOL, C_INT
integer(C_INT32_T), value, intent(in) :: index
logical(C_BOOL), value, intent(in) :: active
integer(C_INT) :: err
end function
function openmc_tally_get_active(index, active) result(err) bind(C)
import C_INT32_T, C_BOOL, C_INT
integer(C_INT32_T), value :: index
logical(C_BOOL), intent(out) :: active
integer(C_INT) :: err
end function
function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C)
import C_INT32_T, C_PTR, C_INT
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: nuclides
integer(C_INT), intent(out) :: n
integer(C_INT) :: err
end function openmc_tally_get_nuclides
function active_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
function active_tallies_data(i) result(tally) bind(C)
import C_INT
integer(C_INT), value :: i
integer(C_INT) :: tally
end function
function active_analog_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
function active_tracklength_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
function active_collision_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
function active_meshsurf_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
function active_surface_tallies_size() result(size) bind(C)
import C_INT
integer(C_INT) :: size
end function
end interface
!===============================================================================
! TALLYOBJECT describes a user-specified tally. The region of phase space to
! tally in is given by the TallyFilters and the results are stored in a
! TallyResult array.
!===============================================================================
type, public :: TallyObject
type(C_PTR) :: ptr
! Basic data
character(len=104) :: name = "" ! user-defined name
real(8) :: volume ! volume of region
integer :: total_score_bins
contains
procedure :: id => tally_get_id
procedure :: set_id => tally_set_id
procedure :: type => tally_get_type
procedure :: set_type => tally_set_type
procedure :: estimator => tally_get_estimator
procedure :: set_estimator => tally_set_estimator
procedure :: depletion_rx => tally_get_depletion_rx
procedure :: n_score_bins => tally_get_n_score_bins
procedure :: score_bins => tally_get_score_bin
procedure :: n_filters => tally_get_n_filters
procedure :: filter => tally_get_filter
procedure :: n_filter_bins => tally_get_n_filter_bins
procedure :: n_nuclide_bins => tally_get_n_nuclide_bins
procedure :: nuclide_bins => tally_get_nuclide_bins
procedure :: energyout_filter => tally_get_energyout_filter
procedure :: deriv => tally_get_deriv
procedure :: set_deriv => tally_set_deriv
end type TallyObject
type, public :: TallyContainer
class(TallyObject), allocatable :: obj
end type TallyContainer
integer(C_INT32_T), public, bind(C) :: n_tallies = 0 ! # of tallies
type(TallyContainer), public, allocatable, target :: tallies(:)
! Dictionary that maps user IDs to indices in 'tallies'
type(DictIntInt), public :: tally_dict
! The largest tally ID that has been specified in the system. This is useful
! in case the code needs to find an ID for a new tally.
integer :: largest_tally_id
contains
function tally_get_id(this) result(t)
class(TallyObject) :: this
integer(C_INT) :: t
interface
function tally_get_id_c(tally) result(t) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: t
end function
end interface
t = tally_get_id_c(this % ptr)
end function
subroutine tally_set_id(this, t)
class(TallyObject) :: this
integer(C_INT) :: t
interface
subroutine tally_set_id_c(tally, t) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: t
end subroutine
end interface
call tally_set_id_c(this % ptr, t)
end subroutine
function tally_get_type(this) result(t)
class(TallyObject) :: this
integer(C_INT) :: t
interface
function tally_get_type_c(tally) result(t) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: t
end function
end interface
t = tally_get_type_c(this % ptr)
end function
subroutine tally_set_type(this, t)
class(TallyObject) :: this
integer(C_INT) :: t
interface
subroutine tally_set_type_c(tally, t) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: t
end subroutine
end interface
call tally_set_type_c(this % ptr, t)
end subroutine
function tally_get_estimator(this) result(e)
class(TallyObject) :: this
integer(C_INT) :: e
interface
function tally_get_estimator_c(tally) result(e) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: e
end function
end interface
e = tally_get_estimator_c(this % ptr)
end function
subroutine tally_set_estimator(this, e)
class(TallyObject) :: this
integer(C_INT) :: e
interface
subroutine tally_set_estimator_c(tally, e) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: e
end subroutine
end interface
call tally_set_estimator_c(this % ptr, e)
end subroutine
function tally_get_depletion_rx(this) result(drx)
class(TallyObject) :: this
logical(C_BOOL) :: drx
interface
function tally_get_depletion_rx_c(tally) result(drx) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: tally
logical(C_BOOl) :: drx
end function
end interface
drx = tally_get_depletion_rx_c(this % ptr)
end function
function tally_get_n_score_bins(this) result(n)
class(TallyObject) :: this
integer(C_INT) :: n
interface
function tally_get_n_scores_c(tally) result(n) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: n
end function
end interface
n = tally_get_n_scores_c(this % ptr)
end function
function tally_get_score_bin(this, i) result(filt)
class(TallyObject) :: this
integer(C_INT) :: i
integer(C_INT32_T) :: filt
interface
function tally_get_score_c(tally, i) result(filt) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: i
integer(C_INT) :: filt
end function
end interface
filt = tally_get_score_c(this % ptr, i-1)
end function
function tally_get_n_filters(this) result(n)
class(TallyObject) :: this
integer(C_INT) :: n
interface
function tally_get_n_filters_c(tally) result(n) bind(C)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), value :: tally
integer(C_INT) :: n
end function
end interface
n = tally_get_n_filters_c(this % ptr)
end function
function tally_get_filter(this, i) result(filt)
class(TallyObject) :: this
integer(C_INT) :: i
integer(C_INT32_T) :: filt
interface
function tally_get_filter_c(tally, i) result(filt) bind(C)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), value :: tally
integer(C_INT), value :: i
integer(C_INT32_T) :: filt
end function
end interface
filt = tally_get_filter_c(this % ptr, i-1)
end function
function tally_get_n_filter_bins(this) result(n_filter_bins)
class(TallyObject) :: this
integer(C_INT32_T) :: n_filter_bins
interface
function tally_get_n_filter_bins_c(tally) result(n_filter_bins) bind(C)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), value :: tally
integer(C_INT32_T) :: n_filter_bins
end function
end interface
n_filter_bins = tally_get_n_filter_bins_c(this % ptr)
end function
function tally_get_n_nuclide_bins(this) result(n)
class(TallyObject) :: this
integer(C_INT) :: n
interface
function tally_get_n_nuclide_bins_c(tally) result(n) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: n
end function
end interface
n = tally_get_n_nuclide_bins_c(this % ptr)
end function
function tally_get_nuclide_bins(this, i) result(nuclide)
class(TallyObject) :: this
integer(C_INT) :: i, nuclide
interface
function tally_get_nuclide_bins_c(tally, i) result(nuclide) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: i
integer(C_INT) :: nuclide
end function
end interface
nuclide = tally_get_nuclide_bins_c(this % ptr, i)
end function
function tally_get_energyout_filter(this) result(filt)
class(TallyObject) :: this
integer(C_INT) :: filt
interface
function tally_get_energyout_filter_c(tally) result(filt) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: filt
end function
end interface
filt = tally_get_energyout_filter_c(this % ptr)
end function
function tally_get_deriv(this) result(deriv)
class(TallyObject) :: this
integer(C_INT) :: deriv
interface
function tally_get_deriv_c(tally) result(deriv) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT) :: deriv
end function
end interface
deriv = tally_get_deriv_c(this % ptr)
end function
subroutine tally_set_deriv(this, deriv)
class(TallyObject) :: this
integer(C_INT) :: deriv
interface
subroutine tally_set_deriv_c(tally, deriv) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally
integer(C_INT), value :: deriv
end subroutine
end interface
call tally_set_deriv_c(this % ptr, deriv)
end subroutine
!===============================================================================
! FREE_MEMORY_TALLY deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_tally()
interface
subroutine free_memory_tally_c() bind(C)
end subroutine free_memory_tally_c
end interface
call free_memory_tally_c()
n_tallies = 0
if (allocated(tallies)) deallocate(tallies)
call tally_dict % clear()
largest_tally_id = 0
end subroutine free_memory_tally
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_tallies(n, index_start, index_end) result(err) bind(C)
! Extend the tallies 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_INT) :: err
integer :: i
type(TallyContainer), allocatable :: temp(:) ! temporary tallies array
interface
subroutine extend_tallies_c() bind(C)
end subroutine
end interface
! Extend the C++ tallies array first
call extend_tallies_c()
if (n_tallies == 0) then
! Allocate tallies array
allocate(tallies(n))
else
! Allocate tallies array with increased size
allocate(temp(n_tallies + n))
! Move original tallies to temporary array
do i = 1, n_tallies
call move_alloc(tallies(i) % obj, temp(i) % obj)
end do
! Move allocation from temporary array
call move_alloc(FROM=temp, TO=tallies)
end if
! Return indices in tallies array
if (present(index_start)) index_start = n_tallies + 1
if (present(index_end)) index_end = n_tallies + n
n_tallies = n_tallies + n
err = 0
end function openmc_extend_tallies
function openmc_get_tally_index(id, index) result(err) bind(C)
! Returns the index in the tallies array of a tally with a given ID
integer(C_INT32_T), value :: id
integer(C_INT32_T), intent(out) :: index
integer(C_INT) :: err
if (allocated(tallies)) then
if (tally_dict % has(id)) then
index = tally_dict % get(id)
err = 0
else
err = E_INVALID_ID
call set_errmsg("No tally exists with ID=" // trim(to_str(id)) // ".")
end if
else
err = E_ALLOCATE
call set_errmsg("Memory has not been allocated for tallies.")
end if
end function openmc_get_tally_index
function openmc_tally_get_estimator(index, estimator) result(err) bind(C)
! Return the type of estimator of a tally
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: estimator
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
estimator = tallies(index) % obj % estimator()
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_get_estimator
function openmc_tally_get_id(index, id) result(err) bind(C)
! Return the ID of a tally
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: id
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
id = tallies(index) % obj % id()
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_get_id
function openmc_tally_set_estimator(index, estimator) result(err) bind(C)
! Set the type of estimator a tally
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR), intent(in) :: estimator(*)
integer(C_INT) :: err
character(:), allocatable :: estimator_
! Convert C string to Fortran string
estimator_ = to_f_string(estimator)
err = 0
if (index >= 1 .and. index <= size(tallies)) then
select case (estimator_)
case ('analog')
call tallies(index) % obj % set_estimator(ESTIMATOR_ANALOG)
case ('tracklength')
call tallies(index) % obj % set_estimator(ESTIMATOR_TRACKLENGTH)
case ('collision')
call tallies(index) % obj % set_estimator(ESTIMATOR_COLLISION)
case default
err = E_INVALID_ARGUMENT
call set_errmsg("Unknown tally estimator: " // trim(estimator_))
end select
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in tally array is out of bounds.")
end if
end function openmc_tally_set_estimator
function openmc_tally_set_id(index, id) result(err) bind(C)
! Set the ID of a tally
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_tallies) then
if (allocated(tallies(index) % obj)) then
if (tally_dict % has(id)) then
call set_errmsg("Two or more tallies use the same unique ID: " &
// to_str(id))
err = E_INVALID_ID
else
call tallies(index) % obj % set_id(id)
call tally_dict % set(id, index)
if (id > largest_tally_id) largest_tally_id = id
err = 0
end if
else
err = E_ALLOCATE
call set_errmsg("Tally type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_set_id
subroutine openmc_get_tally_next_id(id) bind(C)
! Returns an ID number that has not been used by any other tallies.
integer(C_INT32_T), intent(out) :: id
id = largest_tally_id + 1
end subroutine openmc_get_tally_next_id
end module tally_header

12
src/thermal.cpp
View file

@ -587,6 +587,12 @@ ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
}
void free_memory_thermal()
{
data::thermal_scatt.clear();
data::thermal_scatt_map.clear();
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -603,12 +609,6 @@ sab_from_hdf5(hid_t group, const double* temperature, int n)
return data::thermal_scatt.back().get();
}
extern "C" void sab_clear()
{
data::thermal_scatt.clear();
data::thermal_scatt_map.clear();
}
extern "C" bool sab_has_nuclide(int i_sab, const char* name)
{
return data::thermal_scatt[i_sab - 1]->has_nuclide(name);

13
src/volume_calc.cpp
View file

@ -33,7 +33,7 @@ namespace openmc {
//==============================================================================
namespace model {
std::vector<VolumeCalculation> volume_calcs;
std::vector<VolumeCalculation> volume_calcs;
}
//==============================================================================
@ -362,6 +362,11 @@ void VolumeCalculation::check_hit(int i_material, std::vector<int>& indices,
}
}
void free_memory_volume()
{
openmc::model::volume_calcs.clear();
}
} // namespace openmc
//==============================================================================
@ -422,9 +427,3 @@ int openmc_calculate_volumes() {
return 0;
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_volume() { openmc::model::volume_calcs.clear(); }