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

Move CMFD reweight to cpp, minor discrepancies in regression tests

Browse source
This commit is contained in:
Shikhar Kumar 2020-11-24 22:56:31 -05:00
parent 2115cae40a
commit e87e5ba3dc
4 changed files with 240 additions and 166 deletions
Download patch file Download diff file Expand all files Collapse all files

13
include/openmc/capi.h
View file

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

161
openmc/cmfd.py
View file

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

7
openmc/lib/core.py
View file

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

225
src/cmfd_solver.cpp
View file

@ -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<int, 2> indexmap;
int use_all_threads;
RegularMesh* mesh;
std::vector<double> 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<double, 1> count_bank_sites(xt::xtensor<int, 1>& bins, bool* outside)
{
// Determine shape of array for counts
std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng;
std::vector<std::size_t> cnt_shape = {cnt_size};
// Create array of zeros
xt::xarray<double> 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<double>{}.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<double> 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<int, 1> bank_bins({bank_size}, 0);
bool sites_outside;
xt::xtensor<double, 1> sourcecounts = count_bank_sites(bank_bins,
&sites_outside);
// Compute CMFD weightfactors
xt::xtensor<double, 1> weightfactors = xt::xtensor<double, 1>({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<EnergyFilter*>(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<size_t>(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