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

Convert count_bank_sites, Shannon entropy, UFS to C++ (incomplete)

Browse source
This commit is contained in:
Paul Romano 2018-08-30 06:44:37 -05:00
parent fb22413e8d
commit 6a1d653547
22 changed files with 438 additions and 325 deletions
Download patch file Download diff file Expand all files Collapse all files

1
include/openmc/capi.h
View file

@ -37,6 +37,7 @@ extern "C" {
int openmc_filter_set_type(int32_t index, const char* type);
int openmc_finalize();
int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
int openmc_fission_bank(struct Bank** ptr, int64_t* n);
int openmc_get_cell_index(int32_t id, int32_t* index);
int openmc_get_filter_index(int32_t id, int32_t* index);
void openmc_get_filter_next_id(int32_t* id);

37
include/openmc/eigenvalue.h Normal file
View file

@ -0,0 +1,37 @@
#ifndef OPENMC_EIGENVALUE_H
#define OPENMC_EIGENVALUE_H
#include <cstdint> // for int64_t
#include "openmc/particle.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
extern std::vector<double> entropy; //!< Shannon entropy at each generation
extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
extern "C" int64_t n_bank;
//==============================================================================
// Non-member functions
//==============================================================================
//! Calculates the Shannon entropy of the fission source distribution to assess
//! source convergence
extern "C" void shannon_entropy();
//! Determines the source fraction in each UFS mesh cell and reweights the
//! source bank so that the sum of the weights is equal to n_particles. The
//! 'source_frac' variable is used later to bias the production of fission sites
extern "C" void ufs_count_sites();
//! Get UFS weight corresponding to particle's location
extern "C" double ufs_get_weight(const Particle* p);
} // namespace openmc
#endif // OPENMC_EIGENVALUE_H

3
include/openmc/mesh.h
View file

@ -34,6 +34,9 @@ public:
bool intersects(Position r0, Position r1);
void to_hdf5(hid_t group);
xt::xarray<double> count_sites(int64_t n, const Bank* bank,
int n_energy, const double* energies, bool* outside);
int id_ {-1}; //!< User-specified ID
int n_dimension_;
double volume_frac_;

1
include/openmc/message_passing.h
View file

@ -10,6 +10,7 @@ namespace mpi {
extern int rank;
extern int n_procs;
extern bool master;
#ifdef OPENMC_MPI
extern MPI_Datatype bank;

4
src/api.F90
View file

@ -11,7 +11,6 @@ module openmc_api
use hdf5_interface
use material_header
use math
use mesh_header, only: free_memory_mesh
use message_passing
use nuclide_header
use initialize, only: openmc_init_f
@ -305,6 +304,9 @@ contains
interface
subroutine free_memory_source() bind(C)
end subroutine
subroutine free_memory_mesh() bind(C)
end subroutine free_memory_mesh
end interface
call free_memory_geometry()

18
src/bank_header.F90
View file

@ -28,7 +28,7 @@ module bank_header
type(Bank), allocatable, target :: master_fission_bank(:)
#endif
integer(8) :: n_bank ! # of sites in fission bank
integer(C_INT64_T), bind(C) :: n_bank ! # of sites in fission bank
!$omp threadprivate(fission_bank, n_bank)
@ -71,4 +71,20 @@ contains
end if
end function openmc_source_bank
function openmc_fission_bank(ptr, n) result(err) bind(C)
! Return a pointer to the source bank
type(C_PTR), intent(out) :: ptr
integer(C_INT64_T), intent(out) :: n
integer(C_INT) :: err
if (.not. allocated(fission_bank)) then
err = E_ALLOCATE
call set_errmsg("Fission bank has not been allocated.")
else
err = 0
ptr = C_LOC(fission_bank)
n = size(fission_bank)
end if
end function openmc_fission_bank
end module bank_header

95
src/eigenvalue.F90
View file

@ -294,46 +294,6 @@ contains
end subroutine synchronize_bank
!===============================================================================
! SHANNON_ENTROPY calculates the Shannon entropy of the fission source
! distribution to assess source convergence
!===============================================================================
subroutine shannon_entropy()
integer :: i ! index for mesh elements
real(8) :: entropy_gen ! entropy at this generation
logical :: sites_outside ! were there sites outside entropy box?
associate (m => meshes(index_entropy_mesh))
! count number of fission sites over mesh
call count_bank_sites(m, fission_bank, entropy_p, &
size_bank=n_bank, sites_outside=sites_outside)
! display warning message if there were sites outside entropy box
if (sites_outside) then
if (master) call warning("Fission source site(s) outside of entropy box.")
end if
! sum values to obtain shannon entropy
if (master) then
! Normalize to total weight of bank sites
entropy_p = entropy_p / sum(entropy_p)
entropy_gen = ZERO
do i = 1, size(entropy_p, 2)
if (entropy_p(1,i) > ZERO) then
entropy_gen = entropy_gen - &
entropy_p(1,i) * log(entropy_p(1,i))/log(TWO)
end if
end do
! Add value to vector
call entropy % push_back(entropy_gen)
end if
end associate
end subroutine shannon_entropy
!===============================================================================
! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto
! the master processor and normalizes them. This should work whether or not the
@ -577,61 +537,6 @@ contains
end function openmc_get_keff
!===============================================================================
! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and
! reweights the source bank so that the sum of the weights is equal to
! n_particles. The 'source_frac' variable is used later to bias the production
! of fission sites
!===============================================================================
subroutine count_source_for_ufs()
real(8) :: total ! total weight in source bank
logical :: sites_outside ! were there sites outside the ufs mesh?
#ifdef OPENMC_MPI
integer :: n ! total number of ufs mesh cells
integer :: mpi_err ! MPI error code
#endif
associate (m => meshes(index_ufs_mesh))
if (current_batch == 1 .and. current_gen == 1) then
! On the first generation, just assume that the source is already evenly
! distributed so that effectively the production of fission sites is not
! biased
source_frac = m % volume_frac
else
! count number of source sites in each ufs mesh cell
call count_bank_sites(m, source_bank, source_frac, &
sites_outside=sites_outside, size_bank=work)
! Check for sites outside of the mesh
if (master .and. sites_outside) then
call fatal_error("Source sites outside of the UFS mesh!")
end if
#ifdef OPENMC_MPI
! Send source fraction to all processors
n = product(m % dimension)
call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
#endif
! Normalize to total weight to get fraction of source in each cell
total = sum(source_frac)
source_frac = source_frac / total
! Since the total starting weight is not equal to n_particles, we need to
! renormalize the weight of the source sites
source_bank % wgt = source_bank % wgt * n_particles / total
end if
end associate
end subroutine count_source_for_ufs
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission

151
src/eigenvalue.cpp Normal file
View file

@ -0,0 +1,151 @@
#include "openmc/eigenvalue.h"
#include "xtensor/xmath.hpp"
#include "xtensor/xtensor.hpp"
#include "xtensor/xview.hpp"
#include "openmc/capi.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
std::vector<double> entropy;
xt::xtensor<double, 1> source_frac;
//==============================================================================
// Non-member functions
//==============================================================================
void shannon_entropy()
{
// Get pointer to entropy mesh
auto& m = meshes[settings::index_entropy_mesh];
// Get pointer to fission bank
Bank* fission_bank;
int64_t n;
openmc_fission_bank(&fission_bank, &n);
// Get source weight in each mesh bin
bool sites_outside;
xt::xtensor<double, 1> p = m->count_sites(
n_bank, fission_bank, 0, nullptr, &sites_outside);
// display warning message if there were sites outside entropy box
if (sites_outside) {
if (mpi::master) warning("Fission source site(s) outside of entropy box.");
}
// sum values to obtain shannon entropy
if (mpi::master) {
// Normalize to total weight of bank sites
p /= xt::sum(p);
double H = 0.0;
for (auto p_i : p) {
if (p_i > 0.0) {
H -= p_i * std::log(p_i)/std::log(2.0);
}
}
// Add value to vector
entropy.push_back(H);
}
}
void ufs_count_sites()
{
auto &m = meshes[settings::index_ufs_mesh];
if (openmc_current_batch == 1 && openmc_current_gen == 1) {
// On the first generation, just assume that the source is already evenly
// distributed so that effectively the production of fission sites is not
// biased
auto s = xt::view(source_frac, xt::all());
s = m->volume_frac_;
} else {
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
// count number of source sites in each ufs mesh cell
bool sites_outside;
source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr,
&sites_outside);
// Check for sites outside of the mesh
if (mpi::master && sites_outside) {
fatal_error("Source sites outside of the UFS mesh!");
}
#ifdef OPENMC_MPI
// Send source fraction to all processors
int n = xt::prod(m->shape_)();
MPI_Bcast(source_frac.data(), n, MPI_DOUBLE, 0, mpi::intracomm)
#endif
// Normalize to total weight to get fraction of source in each cell
double total = xt::sum(source_frac)();
source_frac /= total;
// Since the total starting weight is not equal to n_particles, we need to
// renormalize the weight of the source sites
for (int i = 0; i < openmc_work; ++i) {
source_bank[i].wgt *= settings::n_particles / total;
}
}
}
double ufs_get_weight(const Particle* p)
{
auto& m = meshes[settings::index_entropy_mesh];
// Determine indices on ufs mesh for current location
// TODO: off by one
int mesh_bin = m->get_bin({p->coord[0].xyz}) - 1;
if (mesh_bin < 0) {
p->write_restart();
fatal_error("Source site outside UFS mesh!");
}
if (source_frac(mesh_bin) != 0.0) {
return m->volume_frac_ / source_frac(mesh_bin);
} else {
return 1.0;
}
}
extern "C" void entropy_to_hdf5(hid_t group)
{
if (settings::entropy_on) {
write_dataset(group, "entropy", entropy);
}
}
extern "C" void entropy_from_hdf5(hid_t group)
{
if (settings::entropy_on) {
read_dataset(group, "entropy", entropy);
}
}
extern "C" double entropy_c(int i)
{
return entropy.at(i - 1);
}
} // namespace openmc

16
src/initialize.cpp
View file

@ -61,7 +61,7 @@ namespace openmc {
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm)
{
openmc::mpi::intracomm = intracomm;
mpi::intracomm = intracomm;
// Initialize MPI
int flag;
@ -69,13 +69,13 @@ void initialize_mpi(MPI_Comm intracomm)
if (!flag) MPI_Init(nullptr, nullptr);
// Determine number of processes and rank for each
MPI_Comm_size(intracomm, &openmc::mpi::n_procs);
MPI_Comm_rank(intracomm, &openmc::mpi::rank);
MPI_Comm_size(intracomm, &mpi::n_procs);
MPI_Comm_rank(intracomm, &mpi::rank);
// Set variable for Fortran side
openmc_n_procs = openmc::mpi::n_procs;
openmc_rank = openmc::mpi::rank;
openmc_master = (openmc::mpi::rank == 0);
openmc_n_procs = mpi::n_procs;
openmc_rank = mpi::rank;
openmc_master = mpi::master = (mpi::rank == 0);
// Create bank datatype
Bank b;
@ -88,8 +88,8 @@ void initialize_mpi(MPI_Comm intracomm)
};
int blocks[] {1, 3, 3, 1, 1};
MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
MPI_Type_create_struct(5, blocks, disp, types, &openmc::mpi::bank);
MPI_Type_commit(&openmc::mpi::bank);
MPI_Type_create_struct(5, blocks, disp, types, &mpi::bank);
MPI_Type_commit(&mpi::bank);
}
#endif // OPENMC_MPI

104
src/input_xml.F90
View file

@ -249,110 +249,6 @@ contains
track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
end if
! Read meshes
call get_node_list(root, "mesh", node_mesh_list)
! Check for user meshes and allocate
n = size(node_mesh_list)
if (n > 0) then
err = openmc_extend_meshes(n, i_start, i_end)
end if
do i = 1, n
associate (m => meshes(i_start + i - 1))
! Instantiate mesh from XML node
call m % from_xml(node_mesh_list(i))
! Add mesh to dictionary
call mesh_dict % set(m % id, i_start + i - 1)
end associate
end do
! Shannon Entropy mesh
if (check_for_node(root, "entropy_mesh")) then
call get_node_value(root, "entropy_mesh", temp_int)
if (mesh_dict % has(temp_int)) then
index_entropy_mesh = mesh_dict % get(temp_int)
else
call fatal_error("Mesh " // to_str(temp_int) // " specified for &
&Shannon entropy does not exist.")
end if
elseif (check_for_node(root, "entropy")) then
call warning("Specifying a Shannon entropy mesh via the <entropy> element &
&is deprecated. Please create a mesh using <mesh> and then reference &
&it by specifying its ID in an <entropy_mesh> element.")
! Get pointer to entropy node
node_entropy = root % child("entropy")
err = openmc_extend_meshes(1, index_entropy_mesh)
associate (m => meshes(index_entropy_mesh))
! Assign ID
m % id = 10000
call m % from_xml(node_entropy)
end associate
end if
if (index_entropy_mesh > 0) then
associate(m => meshes(index_entropy_mesh))
if (.not. allocated(m % dimension)) then
! If the user did not specify how many mesh cells are to be used in
! each direction, we automatically determine an appropriate number of
! cells
m % n_dimension = 3
allocate(m % dimension(3))
m % dimension = ceiling((n_particles/20)**(ONE/THREE))
! Calculate width
m % width = (m % upper_right - m % lower_left) / m % dimension
end if
! Allocate space for storing number of fission sites in each mesh cell
allocate(entropy_p(1, product(m % dimension)))
end associate
! Turn on Shannon entropy calculation
entropy_on = .true.
end if
! Uniform fission source weighting mesh
if (check_for_node(root, "ufs_mesh")) then
call get_node_value(root, "ufs_mesh", temp_int)
if (mesh_dict % has(temp_int)) then
index_ufs_mesh = mesh_dict % get(temp_int)
else
call fatal_error("Mesh " // to_str(temp_int) // " specified for &
&uniform fission site method does not exist.")
end if
elseif (check_for_node(root, "uniform_fs")) then
call warning("Specifying a UFS mesh via the <uniform_fs> element &
&is deprecated. Please create a mesh using <mesh> and then reference &
&it by specifying its ID in a <ufs_mesh> element.")
! Get pointer to ufs node
node_ufs = root % child("uniform_fs")
err = openmc_extend_meshes(1, index_ufs_mesh)
! Allocate mesh object and coordinates on mesh
associate (m => meshes(index_ufs_mesh))
! Assign ID
m % id = 10001
call m % from_xml(node_ufs)
end associate
end if
if (index_ufs_mesh > 0) then
! Allocate array to store source fraction for UFS
allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension)))
! Turn on uniform fission source weighting
ufs = .true.
end if
! Check if the user has specified to write state points
if (check_for_node(root, "state_point")) then

123
src/mesh.cpp
View file

@ -1,9 +1,14 @@
#include "openmc/mesh.h"
#include <algorithm> // for copy
#include <cstddef> // for size_t
#include <cmath> // for ceil
#include <string>
#ifdef OPENMC_MPI
#include "mpi.h"
#endif
#include "xtensor/xbuilder.hpp"
#include "xtensor/xeval.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xsort.hpp"
@ -13,6 +18,8 @@
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/search.h"
#include "openmc/xml_interface.h"
namespace openmc {
@ -521,6 +528,76 @@ void RegularMesh::to_hdf5(hid_t group)
close_group(mesh_group);
}
xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
int n_energy, const double* energies, bool* outside)
{
// Determine shape of array for counts
int m = xt::prod(shape_)();
std::vector<std::size_t> shape;
if (n_energy > 0) {
shape = {m, n_energy};
} else {
shape = {m};
}
// Create array of zeros
xt::xarray<double> cnt {shape, 0.0};
bool outside_ = false;
for (int64_t i = 0; i < n; ++i) {
// determine scoring bin for entropy mesh
int mesh_bin = get_bin({bank[i].xyz});
// if outside mesh, skip particle
if (mesh_bin < 0) {
outside_ = true;
continue;
}
if (n_energy > 0) {
double E = bank[i].E;
int e_bin;
if (E >= energies[0] && E <= energies[n_energy - 1]) {
// determine energy bin
int e_bin = lower_bound_index(energies, energies + n_energy, E);
// Add to appropriate bin
cnt(mesh_bin, e_bin) += bank[i].wgt;
}
} else {
// Add to appropriate bin
cnt(mesh_bin) += bank[i].wgt;
}
}
// Create copy of count data
int total = cnt.size();
double* cnt_reduced = new double[total];
#ifdef OPENMC_MPI
// collect values from all processors
MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0,
mpi::intracomm);
if (outside) {
MPI_Reduce(&outside_, outside, 1, MPI_BOOL, MPI_LOR, 0, mpi::intracomm);
}
// Check if there were sites outside the mesh for any processor
MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
mpi_intracomm, mpi_err)
end if
#else
std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
if (outside) *outside = outside_;
#endif
// Adapt reduced values in array back into an xarray
auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
xt::xarray<double> counts = arr;
return counts;
}
//==============================================================================
// C API functions
//==============================================================================
@ -660,4 +737,50 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
return 0;
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" {
int32_t mesh_id(RegularMesh* m) { return m->id_; }
double mesh_volume_frac(RegularMesh* m) { return m->volume_frac_; }
int mesh_n_dimension(RegularMesh* m) { return m->n_dimension_; }
double* mesh_lower_left(RegularMesh* m) { return m->lower_left_.data(); }
int mesh_get_bin(RegularMesh* m, const double* xyz)
{
return m->get_bin({xyz});
}
void meshes_to_hdf5(hid_t group)
{
// Write number of meshes
hid_t meshes_group = create_group(group, "meshes");
int32_t n_meshes = meshes.size();
write_attribute(meshes_group, "n_meshes", n_meshes);
if (n_meshes > 0) {
// Write IDs of meshes
std::vector<int> ids;
for (const auto& m : meshes) {
m->to_hdf5(meshes_group);
ids.push_back(m->id_);
}
write_attribute(meshes_group, "ids", ids);
}
close_group(meshes_group);
}
void free_memory_mesh()
{
meshes.clear();
mesh_map.clear();
}
}
} // namespace openmc

71
src/mesh_header.F90
View file

@ -11,7 +11,6 @@ module mesh_header
implicit none
private
public :: free_memory_mesh
public :: openmc_extend_meshes
public :: openmc_get_mesh_index
public :: openmc_mesh_get_id
@ -27,22 +26,23 @@ module mesh_header
!===============================================================================
type, public :: RegularMesh
integer :: id = -1 ! user-specified id
integer :: type = MESH_REGULAR ! rectangular, hexagonal
type(C_PTR) :: ptr
integer(C_INT) :: n_dimension ! rank of mesh
real(8) :: volume_frac ! volume fraction of each cell
integer(C_INT), allocatable :: dimension(:) ! number of cells in each direction
real(C_DOUBLE), allocatable :: lower_left(:) ! lower-left corner of mesh
real(C_DOUBLE), allocatable :: upper_right(:) ! upper-right corner of mesh
real(C_DOUBLE), allocatable :: width(:) ! width of each mesh cell
contains
procedure :: id => regular_id
procedure :: volume_frac => regular_volume_frac
procedure :: from_xml => regular_from_xml
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
procedure :: intersects => regular_intersects
procedure :: to_hdf5 => regular_to_hdf5
end type RegularMesh
integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes
@ -121,6 +121,36 @@ module mesh_header
contains
function regular_id(this) result(id)
class(RegularMesh), intent(in) :: this
integer(C_INT32_T) :: id
interface
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
end interface
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
interface
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
end interface
volume_frac = mesh_volume_frac(this % ptr)
end function
subroutine regular_from_xml(this, node)
class(RegularMesh), intent(inout) :: this
type(XMLNode), intent(in) :: node
@ -598,35 +628,4 @@ contains
end function mesh_intersects_3d
!===============================================================================
! TO_HDF5 writes the mesh data to an HDF5 group
!===============================================================================
subroutine regular_to_hdf5(this, group)
class(RegularMesh), intent(in) :: this
integer(HID_T), intent(in) :: group
integer(HID_T) :: mesh_group
mesh_group = create_group(group, "mesh " // trim(to_str(this % id)))
call write_dataset(mesh_group, "type", "regular")
call write_dataset(mesh_group, "dimension", this % dimension)
call write_dataset(mesh_group, "lower_left", this % lower_left)
call write_dataset(mesh_group, "upper_right", this % upper_right)
call write_dataset(mesh_group, "width", this % width)
call close_group(mesh_group)
end subroutine regular_to_hdf5
!===============================================================================
! FREE_MEMORY_MESH deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_mesh()
n_meshes = 0
if (allocated(meshes)) deallocate(meshes)
call mesh_dict % clear()
end subroutine free_memory_mesh
end module mesh_header

1
src/message_passing.cpp
View file

@ -5,6 +5,7 @@ namespace mpi {
int rank {0};
int n_procs {1};
bool master {false};
#ifdef OPENMC_MPI
MPI_Comm intracomm;

12
src/output.F90
View file

@ -33,6 +33,14 @@ module output
integer :: ou = OUTPUT_UNIT
integer :: eu = ERROR_UNIT
interface
function entropy(i) result(h) bind(C, name='entropy_c')
import C_INT, C_DOUBLE
integer(C_INT), value :: i
real(C_DOUBLE) :: h
end function
end interface
contains
!===============================================================================
@ -335,7 +343,7 @@ contains
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy % data(i)
entropy(i)
if (n > 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
@ -369,7 +377,7 @@ contains
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy % data(i)
entropy(i)
! write out accumulated k-effective if after first active batch
if (n > 1) then

23
src/physics.F90
View file

@ -1187,6 +1187,14 @@ contains
real(8) :: weight ! weight adjustment for ufs method
type(Nuclide), pointer :: nuc
interface
function ufs_get_weight(p) result(weight) bind(C)
import Particle, C_DOUBLE
type(Particle), intent(in) :: p
real(C_DOUBLE) :: WEIGHT
end function
end interface
! Get pointers
nuc => nuclides(i_nuclide)
@ -1196,20 +1204,7 @@ contains
! the expected number of fission sites produced
if (ufs) then
associate (m => meshes(index_ufs_mesh))
! Determine indices on ufs mesh for current location
call m % get_bin(p % coord(1) % xyz, mesh_bin)
if (mesh_bin == NO_BIN_FOUND) then
call particle_write_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (source_frac(1, mesh_bin) /= ZERO) then
weight = m % volume_frac / source_frac(1, mesh_bin)
else
weight = ONE
end if
end associate
weight = ufs_get_weight(p)
else
weight = ONE
end if

24
src/physics_mg.F90
View file

@ -174,27 +174,21 @@ contains
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
interface
function ufs_get_weight(p) result(weight) bind(C)
import Particle, C_DOUBLE
type(Particle), intent(in) :: p
real(C_DOUBLE) :: WEIGHT
end function
end interface
! TODO: Heat generation from fission
! If uniform fission source weighting is turned on, we increase of decrease
! the expected number of fission sites produced
if (ufs) then
associate (m => meshes(index_ufs_mesh))
! Determine indices on ufs mesh for current location
call m % get_bin(p % coord(1) % xyz, mesh_bin)
if (mesh_bin == NO_BIN_FOUND) then
call particle_write_restart(p)
call fatal_error("Source site outside UFS mesh!")
end if
if (source_frac(1, mesh_bin) /= ZERO) then
weight = m % volume_frac / source_frac(1, mesh_bin)
else
weight = ONE
end if
end associate
weight = ufs_get_weight(p)
else
weight = ONE
end if

8
src/settings.cpp
View file

@ -537,10 +537,6 @@ read_settings_xml()
m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_;
}
// TODO: Allocate entropy_P
// Allocate space for storing number of fission sites in each mesh cell
//allocate(entropy_p(1, product(m % dimension)))
// Turn on Shannon entropy calculation
settings::entropy_on = true;
}
@ -574,10 +570,6 @@ read_settings_xml()
}
if (index_ufs_mesh >= 0) {
// Allocate array to store source fraction for UFS
// TODO: Allocate source_frac
//allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension)))
// Turn on uniform fission source weighting
settings::ufs_on = true;
}

12
src/simulation.F90
View file

@ -10,7 +10,7 @@ module simulation
use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd
use cmfd_header, only: cmfd_on
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
use eigenvalue, only: calculate_average_keff, &
calculate_generation_keff, shannon_entropy, &
synchronize_bank, keff_generation, k_sum
#ifdef _OPENMP
@ -234,12 +234,17 @@ contains
subroutine initialize_generation()
interface
subroutine ufs_count_sites() bind(C)
end subroutine
end interface
if (run_mode == MODE_EIGENVALUE) then
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
if (ufs) call ufs_count_sites()
! Store current value of tracklength k
keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH)
@ -256,6 +261,9 @@ contains
interface
subroutine fill_source_bank_fixedsource() bind(C)
end subroutine
subroutine shannon_entropy() bind(C)
end subroutine
end interface
! Update global tallies with the omp private accumulation variables

9
src/simulation_header.F90
View file

@ -47,13 +47,6 @@ module simulation_header
real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
! Shannon entropy
type(VectorReal) :: entropy ! shannon entropy at each generation
real(8), allocatable :: entropy_p(:,:) ! % of source sites in each cell
! Uniform fission source weighting
real(8), allocatable :: source_frac(:,:)
! ============================================================================
! PARALLEL PROCESSING VARIABLES
@ -87,8 +80,6 @@ contains
!===============================================================================
subroutine free_memory_simulation()
if (allocated(entropy_p)) deallocate(entropy_p)
if (allocated(source_frac)) deallocate(source_frac)
if (allocated(work_index)) deallocate(work_index)
call k_generation % clear()

48
src/state_point.F90
View file

@ -20,7 +20,6 @@ module state_point
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
use hdf5_interface
use mesh_header, only: RegularMesh, meshes, n_meshes
use message_passing
use mgxs_interface
use nuclide_header, only: nuclides
@ -79,6 +78,17 @@ contains
character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
logical :: parallel
interface
subroutine meshes_to_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine entropy_to_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
end interface
err = 0
! Set the filename
@ -163,9 +173,7 @@ contains
call write_dataset(file_id, "generations_per_batch", gen_per_batch)
k = k_generation % size()
call write_dataset(file_id, "k_generation", k_generation % data(1:k))
if (entropy_on) then
call write_dataset(file_id, "entropy", entropy % data(1:k))
end if
call entropy_to_hdf5()
call write_dataset(file_id, "k_col_abs", k_col_abs)
call write_dataset(file_id, "k_col_tra", k_col_tra)
call write_dataset(file_id, "k_abs_tra", k_abs_tra)
@ -192,26 +200,8 @@ contains
tallies_group = create_group(file_id, "tallies")
! Write number of meshes
meshes_group = create_group(tallies_group, "meshes")
call write_attribute(meshes_group, "n_meshes", n_meshes)
if (n_meshes > 0) then
! Write IDs of meshes
allocate(id_array(n_meshes))
do i = 1, n_meshes
id_array(i) = meshes(i) % id
end do
call write_attribute(meshes_group, "ids", id_array)
deallocate(id_array)
! Write information for meshes
MESH_LOOP: do i = 1, n_meshes
call meshes(i) % to_hdf5(meshes_group)
end do MESH_LOOP
end if
call close_group(meshes_group)
! Write meshes
call meshes_to_hdf5(tallies_group)
! Write information for derivatives.
if (size(tally_derivs) > 0) then
@ -641,6 +631,11 @@ contains
logical :: source_present
character(MAX_WORD_LEN) :: word
interface
subroutine entropy_from_hdf5() bind(C)
end subroutine
end interface
! Write message
call write_message("Loading state point " // trim(path_state_point) &
// "...", 5)
@ -722,10 +717,7 @@ contains
call k_generation % resize(n)
call read_dataset(k_generation % data(1:n), file_id, "k_generation")
if (entropy_on) then
call entropy % resize(n)
call read_dataset(entropy % data(1:n), file_id, "entropy")
end if
call entropy_from_hdf5()
call read_dataset(k_col_abs, file_id, "k_col_abs")
call read_dataset(k_col_tra, file_id, "k_col_tra")
call read_dataset(k_abs_tra, file_id, "k_abs_tra")

1
src/summary.F90
View file

@ -6,7 +6,6 @@ module summary
use geometry_header
use hdf5_interface
use material_header, only: Material, n_materials, openmc_material_get_volume
use mesh_header, only: RegularMesh
use message_passing
use mgxs_interface
use nuclide_header

1
src/tallies/tally.F90
View file

@ -8,7 +8,6 @@ module tally
use error, only: fatal_error
use geometry_header
use math, only: t_percentile
use mesh_header, only: RegularMesh, meshes
use message_passing
use mgxs_interface
use nuclide_header