From e87e5ba3dcd2ff8f965dc1c2bea39e45e07e9701 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 24 Nov 2020 22:56:31 -0500 Subject: [PATCH] Move CMFD reweight to cpp, minor discrepancies in regression tests --- include/openmc/capi.h | 13 ++- openmc/cmfd.py | 161 +++--------------------------- openmc/lib/core.py | 7 +- src/cmfd_solver.cpp | 225 +++++++++++++++++++++++++++++++++++++++--- 4 files changed, 240 insertions(+), 166 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 6d82929d1..19f31479a 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -131,6 +131,18 @@ extern "C" { int openmc_zernike_filter_set_params(int32_t index, const double* x, const double* y, const double* r); + //! Sets the mesh and energy grid for CMFD reweight + //! \param[in] id of CMFD Mesh Tally + //! \param[in] CMFD normalization factor + //! \param[in] CMFD weight clipping factor + extern "C" void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices, + const double norm); + + //! Sets the mesh and energy grid for CMFD reweight + //! \param[in] whether or not to run CMFD feedback + //! \param[in] computed CMFD source + extern "C" void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src); + //! Sets the fixed variables that are used for CMFD linear solver //! \param[in] CSR format index pointer array of loss matrix //! \param[in] length of indptr @@ -143,7 +155,6 @@ extern "C" { extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr, const int* indices, int n_elements, int dim, double spectral, - const int* cmfd_indices, const int* map, bool use_all_threads); //! Runs a Gauss Seidel linear solver to solve CMFD matrix equations diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 78901b03b..a99b60619 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -366,8 +366,6 @@ class CMFDRun: self._current = None self._cmfd_src = None self._openmc_src = None - self._sourcecounts = None - self._weightfactors = None self._entropy = [] self._balance = [] self._src_cmp = [] @@ -914,7 +912,7 @@ class CMFDRun: args = temp_loss.indptr, len(temp_loss.indptr), \ temp_loss.indices, len(temp_loss.indices), n, \ - self._spectral, self._indices, coremap, self._use_all_threads + self._spectral, coremap, self._use_all_threads return openmc.lib._dll.openmc_initialize_linsolver(*args) def _write_cmfd_output(self): @@ -1171,8 +1169,12 @@ class CMFDRun: # Calculate fission source self._calc_fission_source() - # Calculate weight factors - self._cmfd_reweight() + # Calculate weight factors through C++ and manipulate CMFD + # source into a 1-D vector that matches C++ array ordering + src_flipped = np.flip(self._cmfd_src, axis=3) + src_swapped = np.swapaxes(src_flipped, 0, 2) + args = self._feedback, src_swapped.flatten() + openmc.lib._dll.openmc_cmfd_reweight(*args) # Stop CMFD timer if openmc.lib.master(): @@ -1390,151 +1392,6 @@ class CMFDRun: self._src_cmp.append(np.sqrt(1.0 / self._norm * np.sum((self._cmfd_src - self._openmc_src)**2))) - def _cmfd_reweight(self): - """Performs weighting of particles in source bank""" - # Get spatial dimensions and energy groups - nx, ny, nz, ng = self._indices - - # Count bank site in mesh and reverse due to egrid structured - outside = self._count_bank_sites() - - # Check and raise error if source sites exist outside of CMFD mesh - if openmc.lib.master() and outside: - raise OpenMCError('Source sites outside of the CMFD mesh') - - # Have master compute weight factors, ignore any zeros in - # sourcecounts or cmfd_src - if openmc.lib.master(): - # Compute normalization factor - norm = np.sum(self._sourcecounts) / np.sum(self._cmfd_src) - - # Define target reshape dimensions for sourcecounts. This - # defines how self._sourcecounts is ordered by dimension - target_shape = [nz, ny, nx, ng] - - # Reshape sourcecounts to target shape. Swap x and z axes so - # that the shape is now [nx, ny, nz, ng] - sourcecounts = np.swapaxes( - self._sourcecounts.reshape(target_shape), 0, 2) - - # Flip index of energy dimension - sourcecounts = np.flip(sourcecounts, axis=3) - - # Compute weight factors - div_condition = np.logical_and(sourcecounts > 0, - self._cmfd_src > 0) - self._weightfactors = (np.divide(self._cmfd_src * norm, - sourcecounts, where=div_condition, - out=np.ones_like(self._cmfd_src), - dtype=np.float32)) - - if not self._feedback: - return - - # Broadcast weight factors to all procs - if have_mpi: - self._weightfactors = self._intracomm.bcast( - self._weightfactors) - - m = openmc.lib.meshes[self._mesh_id] - energy = self._egrid - ng = self._indices[3] - - # Get locations and energies of all particles in source bank - source_xyz = openmc.lib.source_bank()['r'] - source_energies = openmc.lib.source_bank()['E'] - - # Convert xyz location to the CMFD mesh index - mesh_ijk = np.floor((source_xyz - m.lower_left)/m.width).astype(int) - - # Determine which energy bin each particle's energy belongs to - # Separate into cases bases on where source energies lies on egrid - energy_bins = np.zeros(len(source_energies), dtype=int) - idx = np.where(source_energies < energy[0]) - energy_bins[idx] = ng - 1 - idx = np.where(source_energies > energy[-1]) - energy_bins[idx] = 0 - idx = np.where((source_energies >= energy[0]) & - (source_energies <= energy[-1])) - energy_bins[idx] = ng - np.digitize(source_energies[idx], energy) - - # Determine weight factor of each particle based on its mesh index - # and energy bin and updates its weight - openmc.lib.source_bank()['wgt'] *= self._weightfactors[ - mesh_ijk[:,0], mesh_ijk[:,1], mesh_ijk[:,2], energy_bins] - - if openmc.lib.master() and np.any(source_energies < energy[0]): - print(' WARNING: Source point below energy grid') - sys.stdout.flush() - if openmc.lib.master() and np.any(source_energies > energy[-1]): - print(' WARNING: Source point above energy grid') - sys.stdout.flush() - - def _count_bank_sites(self): - """Determines the number of fission bank sites in each cell of a given - mesh and energy group structure. - Returns - ------- - bool - Wheter any source sites outside of CMFD mesh were found - - """ - # Initialize variables - m = openmc.lib.meshes[self._mesh_id] - bank = openmc.lib.source_bank() - energy = self._egrid - sites_outside = np.zeros(1, dtype=bool) - nxnynz = np.prod(self._indices[0:3]) - ng = self._indices[3] - - outside = np.zeros(1, dtype=bool) - self._sourcecounts = np.zeros((nxnynz, ng)) - count = np.zeros(self._sourcecounts.shape) - - # Get location and energy of each particle in source bank - source_xyz = openmc.lib.source_bank()['r'] - source_energies = openmc.lib.source_bank()['E'] - - # Convert xyz location to mesh index and ravel index to scalar - mesh_locations = np.floor((source_xyz - m.lower_left) / m.width) - mesh_bins = mesh_locations[:,2] * m.dimension[1] * m.dimension[0] + \ - mesh_locations[:,1] * m.dimension[0] + mesh_locations[:,0] - - # Check if any source locations lie outside of defined CMFD mesh - if np.any(mesh_bins < 0) or np.any(mesh_bins >= np.prod(m.dimension)): - outside[0] = True - - # Determine which energy bin each particle's energy belongs to - # Separate into cases bases on where source energies lies on egrid - energy_bins = np.zeros(len(source_energies), dtype=int) - idx = np.where(source_energies < energy[0]) - energy_bins[idx] = 0 - idx = np.where(source_energies > energy[-1]) - energy_bins[idx] = ng - 1 - idx = np.where((source_energies >= energy[0]) & - (source_energies <= energy[-1])) - energy_bins[idx] = np.digitize(source_energies[idx], energy) - 1 - - # Determine all unique combinations of mesh bin and energy bin, and - # count number of particles that belong to these combinations - idx, counts = np.unique(np.array([mesh_bins, energy_bins]), axis=1, - return_counts=True) - - # Store counts to appropriate mesh-energy combination - count[idx[0].astype(int), idx[1].astype(int)] = counts - - if have_mpi: - # Collect values of count from all processors - self._intracomm.Reduce(count, self._sourcecounts, MPI.SUM) - # Check if there were sites outside the mesh for any processor - self._intracomm.Reduce(outside, sites_outside, MPI.LOR) - # Deal with case if MPI not defined (only one proc) - else: - sites_outside = outside - self._sourcecounts = count - - return sites_outside[0] - def _build_loss_matrix(self, adjoint): # Extract spatial and energy indices and define matrix dimension ng = self._indices[3] @@ -3045,3 +2902,7 @@ class CMFDRun: # Set all tallies to be active from beginning cmfd_tally.active = True + + # Initialize CMFD mesh and energy grid in C++ for CMFD reweight + args = self._tally_ids[0], self._indices, self._norm + openmc.lib._dll.openmc_initialize_mesh_egrid(*args) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 3fe75dccf..f24f5d7dc 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -31,6 +31,8 @@ _array_1d_dble = np.ctypeslib.ndpointer(dtype=np.double, ndim=1, _dll.openmc_calculate_volumes.restype = c_int _dll.openmc_calculate_volumes.errcheck = _error_handler +_dll.openmc_cmfd_reweight.argtypes = c_bool, _array_1d_dble +_dll.openmc_cmfd_reweight.restype = None _dll.openmc_finalize.restype = c_int _dll.openmc_finalize.errcheck = _error_handler _dll.openmc_find_cell.argtypes = [POINTER(c_double*3), POINTER(c_int32), @@ -45,8 +47,11 @@ _dll.openmc_init.errcheck = _error_handler _dll.openmc_get_keff.argtypes = [POINTER(c_double*2)] _dll.openmc_get_keff.restype = c_int _dll.openmc_get_keff.errcheck = _error_handler +_init_mesh_egrid_args = [c_int, _array_1d_int, c_double] +_dll.openmc_initialize_mesh_egrid.argtypes = _init_mesh_egrid_args +_dll.openmc_initialize_mesh_egrid.restype = None _init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int, - c_double, _array_1d_int, _array_1d_int, c_bool] + c_double, _array_1d_int, c_bool] _dll.openmc_initialize_linsolver.argtypes = _init_linsolver_argtypes _dll.openmc_initialize_linsolver.restype = None _dll.openmc_is_statepoint_batch.restype = c_bool diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 8406c9899..a3f354c17 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -8,9 +8,15 @@ #endif #include "xtensor/xtensor.hpp" +#include "openmc/bank.h" #include "openmc/error.h" #include "openmc/constants.h" #include "openmc/capi.h" +#include "openmc/mesh.h" +#include "openmc/message_passing.h" +#include "openmc/tallies/filter_energy.h" +#include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/tally.h" namespace openmc { @@ -34,8 +40,203 @@ xt::xtensor indexmap; int use_all_threads; +RegularMesh* mesh; + +std::vector egrid; + +double norm; + } // namespace cmfd +//============================================================================== +// GET_CMFD_ENERGY_BIN returns the energy bin for a source site energy +//============================================================================== + +int get_cmfd_energy_bin(const double E) +{ + // Check if energy is out of grid bounds + if (E < cmfd::egrid[0]) { + // throw warning message + warning("Detected source point below energy grid"); + return 0; + } else if (E >= cmfd::egrid[cmfd::ng]) { + // throw warning message + warning("Detected source point above energy grid"); + return cmfd::ng - 1; + } else { + // Iterate through energy grid to find matching bin + for (int g = 0; g < cmfd::ng; g++) { + if (E >= cmfd::egrid[g] && E < cmfd::egrid[g+1]) { + return g; + } + } + } + // Return -1 by default + return -1; +} + +//============================================================================== +// COUNT_BANK_SITES bins fission sites according to CMFD mesh and energy +//============================================================================== + +xt::xtensor count_bank_sites(xt::xtensor& bins, bool* outside) +{ + // Determine shape of array for counts + std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng; + std::vector cnt_shape = {cnt_size}; + + // Create array of zeros + xt::xarray cnt {cnt_shape, 0.0}; + bool outside_ = false; + + auto bank_size = simulation::source_bank.size(); + for (int i = 0; i < bank_size; i++) { + const auto& site = simulation::source_bank[i]; + + // determine scoring bin for CMFD mesh + int mesh_bin = cmfd::mesh->get_bin(site.r); + + // if outside mesh, skip particle + if (mesh_bin < 0) { + outside_ = true; + continue; + } + + // determine scoring bin for CMFD energy + int energy_bin = get_cmfd_energy_bin(site.E); + + // add to appropriate bin + cnt(mesh_bin*cmfd::ng+energy_bin) += site.wgt; + + // store bin index which is used again when updating weights + bins[i] = mesh_bin*cmfd::ng+energy_bin; + } + + // Create copy of count data. Since ownership will be acquired by xtensor, + // std::allocator must be used to avoid Valgrind mismatched free() / delete + // warnings. + int total = cnt.size(); + double* cnt_reduced = std::allocator{}.allocate(total); + +#ifdef OPENMC_MPI + // collect values from all processors + MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, + mpi::intracomm); + + // Check if there were sites outside the mesh for any processor + MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm); + +#else + std::copy(cnt.data(), cnt.data() + total, cnt_reduced); + *outside = outside_; +#endif + + // Adapt reduced values in array back into an xarray + auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), cnt_shape); + xt::xarray counts = arr; + + return counts; +} + +//============================================================================== +// OPENMC_CMFD_REWEIGHT performs reweighting of particles in source bank +//============================================================================== + +extern "C" +void openmc_cmfd_reweight(const bool feedback, const double* cmfd_src) +{ + // Get size of source bank and cmfd_src + auto bank_size = simulation::source_bank.size(); + std::size_t src_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng; + + // count bank sites for CMFD mesh, store bins in bank_bins for reweighting + xt::xtensor bank_bins({bank_size}, 0); + bool sites_outside; + xt::xtensor sourcecounts = count_bank_sites(bank_bins, + &sites_outside); + + // Compute CMFD weightfactors + xt::xtensor weightfactors = xt::xtensor({src_size}, 1.); + if (mpi::master) { + if (sites_outside) { + fatal_error("Source sites outside of the CMFD mesh"); + } + + double norm = xt::sum(sourcecounts)()/cmfd::norm; + for (int i = 0; i < src_size; i++) { + if (sourcecounts[i] > 0 && cmfd_src[i] > 0) { + weightfactors[i] = cmfd_src[i] * norm / sourcecounts[i]; + } + } + } + + if (!feedback) { + return; + } + +#ifdef OPENMC_MPI + // Send weightfactors to all processors + MPI_Bcast(weightfactors.data(), src_size, MPI_DOUBLE, 0, mpi::intracomm); +#endif + + // Iterate through fission bank and update particle weights + for (int64_t i = 0; i < bank_size; i++) { + auto& site = simulation::source_bank[i]; + site.wgt *= weightfactors(bank_bins(i)); + } +} + +//============================================================================== +// OPENMC_INITIALIZE_MESH_EGRID sets the mesh and energy grid for CMFD reweight +//============================================================================== + +extern "C" +void openmc_initialize_mesh_egrid(const int meshtally_id, const int* cmfd_indices, + const double norm) +{ + // Make sure all CMFD memory is freed + free_memory_cmfd(); + + // Set CMFD indices + cmfd::nx = cmfd_indices[0]; + cmfd::ny = cmfd_indices[1]; + cmfd::nz = cmfd_indices[2]; + cmfd::ng = cmfd_indices[3]; + + // Set CMFD reweight properties + cmfd::norm = norm; + + // Find index corresponding to tally id + int32_t tally_index; + openmc_get_tally_index(meshtally_id, &tally_index); + + // Get filters assocaited with tally + auto tally_filters = model::tallies[tally_index]->filters(); + + // Get mesh filter index + auto meshfilter_index = tally_filters[0]; + + // Store energy filter index if defined, otherwise set to -1 + auto energy_index = (tally_filters.size() == 2) ? tally_filters[1] : -1; + + // Get mesh index from mesh filter index + int32_t mesh_index; + openmc_mesh_filter_get_mesh(meshfilter_index, &mesh_index); + + // Get mesh from mesh index + cmfd::mesh = get_regular_mesh(mesh_index); + + // Get energy bins from energy index, otherwise use default + if (energy_index != -1) + { + auto efilt_base = model::tally_filters[energy_index].get(); + auto* efilt = dynamic_cast(efilt_base); + cmfd::egrid = efilt->bins(); + } else { + cmfd::egrid = {0.0, INFTY}; + } +} + //============================================================================== // MATRIX_TO_INDICES converts a matrix index to spatial and group // indices @@ -309,12 +510,9 @@ int cmfd_linsolver_ng(const double* A_data, const double* b, double* x, extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr, const int* indices, int n_elements, int dim, - double spectral, const int* cmfd_indices, - const int* map, bool use_all_threads) + double spectral, const int* map, + bool use_all_threads) { - // Make sure vectors are empty - free_memory_cmfd(); - // Store elements of indptr for (int i = 0; i < len_indptr; i++) cmfd::indptr.push_back(indptr[i]); @@ -327,15 +525,8 @@ void openmc_initialize_linsolver(const int* indptr, int len_indptr, cmfd::dim = dim; cmfd::spectral = spectral; - // Set number of groups - cmfd::ng = cmfd_indices[3]; - - // Set problem dimensions and indexmap if 1 or 2 group problem + // Set indexmap if 1 or 2 group problem if (cmfd::ng == 1 || cmfd::ng == 2) { - cmfd::nx = cmfd_indices[0]; - cmfd::ny = cmfd_indices[1]; - cmfd::nz = cmfd_indices[2]; - // Resize indexmap and set its elements cmfd::indexmap.resize({static_cast(dim), 3}); set_indexmap(map); @@ -366,10 +557,16 @@ int openmc_run_linsolver(const double* A_data, const double* b, double* x, void free_memory_cmfd() { + // Clear std::vectors cmfd::indptr.clear(); cmfd::indices.clear(); - // Resize indexmap to be an empty array + cmfd::egrid.clear(); + + // Resize xtensors to be empty cmfd::indexmap.resize({0}); + + // Set pointers to null + cmfd::mesh = nullptr; } } // namespace openmc