diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 7ce42aa47..c155e7bd1 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -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); diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h new file mode 100644 index 000000000..c2e76e26f --- /dev/null +++ b/include/openmc/eigenvalue.h @@ -0,0 +1,37 @@ +#ifndef OPENMC_EIGENVALUE_H +#define OPENMC_EIGENVALUE_H + +#include // for int64_t + +#include "openmc/particle.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +extern std::vector entropy; //!< Shannon entropy at each generation +extern xt::xtensor 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 diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 71649fa61..51a0a627f 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -34,6 +34,9 @@ public: bool intersects(Position r0, Position r1); void to_hdf5(hid_t group); + xt::xarray 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_; diff --git a/include/openmc/message_passing.h b/include/openmc/message_passing.h index 1dab43139..c910210db 100644 --- a/include/openmc/message_passing.h +++ b/include/openmc/message_passing.h @@ -10,6 +10,7 @@ namespace mpi { extern int rank; extern int n_procs; + extern bool master; #ifdef OPENMC_MPI extern MPI_Datatype bank; diff --git a/src/api.F90 b/src/api.F90 index 2a6775a0c..b39a2f8fa 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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() diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 2a998740e..c39bf8b78 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -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 diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index d47bc4dab..38284c279 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp new file mode 100644 index 000000000..9ce5b9ef8 --- /dev/null +++ b/src/eigenvalue.cpp @@ -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 entropy; +xt::xtensor 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 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 diff --git a/src/initialize.cpp b/src/initialize.cpp index 2c6840151..19fc82c85 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ffee75450..61b320804 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 element & - &is deprecated. Please create a mesh using and then reference & - &it by specifying its ID in an 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 element & - &is deprecated. Please create a mesh using and then reference & - &it by specifying its ID in a 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 diff --git a/src/mesh.cpp b/src/mesh.cpp index 8aecd1a22..47f03ad3d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,9 +1,14 @@ #include "openmc/mesh.h" +#include // for copy #include // for size_t #include // for ceil #include +#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 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 shape; + if (n_energy > 0) { + shape = {m, n_energy}; + } else { + shape = {m}; + } + + // Create array of zeros + xt::xarray 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 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 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 diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 52c220976..8a8d16c08 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -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 diff --git a/src/message_passing.cpp b/src/message_passing.cpp index aac677c6e..ec2fc2d68 100644 --- a/src/message_passing.cpp +++ b/src/message_passing.cpp @@ -5,6 +5,7 @@ namespace mpi { int rank {0}; int n_procs {1}; +bool master {false}; #ifdef OPENMC_MPI MPI_Comm intracomm; diff --git a/src/output.F90 b/src/output.F90 index b7b418a39..e28695667 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 992e48767..f7d80261a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 40df6e1cb..84695153c 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -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 diff --git a/src/settings.cpp b/src/settings.cpp index 86f3c1fef..7b00b9eca 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -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; } diff --git a/src/simulation.F90 b/src/simulation.F90 index d334b46b8..8b7cfdb27 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -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 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 0e38b76f2..9b7c68cd4 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -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() diff --git a/src/state_point.F90 b/src/state_point.F90 index bc6474de2..b7f4e1377 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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") diff --git a/src/summary.F90 b/src/summary.F90 index 750a9899f..0fafd7e5d 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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 diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 9967606d1..a618f351e 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -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