Merge branch 'develop' into cell-instance-filter

This commit is contained in:
Paul Romano 2020-01-09 07:26:43 -06:00
commit 9787e9e8c0
37 changed files with 646 additions and 433 deletions

View file

@ -1,4 +1,4 @@
Copyright (c) 2011-2019 Massachusetts Institute of Technology and OpenMC contributors
Copyright (c) 2011-2020 Massachusetts Institute of Technology and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

View file

@ -72,7 +72,7 @@ master_doc = 'index'
# General information about the project.
project = 'OpenMC'
copyright = '2011-2019, Massachusetts Institute of Technology and OpenMC contributors'
copyright = '2011-2020, Massachusetts Institute of Technology and OpenMC contributors'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the

View file

@ -435,16 +435,19 @@ attributes/sub-elements:
:type:
The type of spatial distribution. Valid options are "box", "fission",
"point", "cartesian", and "spherical". A "box" spatial distribution has
coordinates sampled uniformly in a parallelepiped. A "fission" spatial
distribution samples locations from a "box" distribution but only
locations in fissionable materials are accepted. A "point" spatial
distribution has coordinates specified by a triplet. An "cartesian"
spatial distribution specifies independent distributions of x-, y-, and
z-coordinates. A "spherical" spatial distribution specifies independent
distributions of r-, theta-, and phi-coordinates where theta is the angle
with respect to the z-axis, phi is the azimuthal angle, and the sphere is
centered on the coordinate (x0,y0,z0).
"point", "cartesian", "cylindrical", and "spherical". A "box" spatial
distribution has coordinates sampled uniformly in a parallelepiped. A
"fission" spatial distribution samples locations from a "box"
distribution but only locations in fissionable materials are accepted.
A "point" spatial distribution has coordinates specified by a triplet.
A "cartesian" spatial distribution specifies independent distributions of
x-, y-, and z-coordinates. A "cylindrical" spatial distribution specifies
independent distributions of r-, phi-, and z-coordinates where phi is the
azimuthal angle and the origin for the cylindrical coordinate system is
specified by origin. A "spherical" spatial distribution specifies
independent distributions of r-, theta-, and phi-coordinates where theta
is the angle with respect to the z-axis, phi is the azimuthal angle, and
the sphere is centered on the coordinate (x0,y0,z0).
*Default*: None
@ -462,6 +465,9 @@ attributes/sub-elements:
For an "cartesian" distribution, no parameters are specified. Instead,
the ``x``, ``y``, and ``z`` elements must be specified.
For a "cylindrical" distribution, no parameters are specified. Instead,
the ``r``, ``phi``, ``z``, and ``origin`` elements must be specified.
For a "spherical" distribution, no parameters are specified. Instead,
the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified.
@ -480,15 +486,16 @@ attributes/sub-elements:
:ref:`univariate`).
:z:
For an "cartesian" distribution, this element specifies the distribution
of z-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
:ref:`univariate`).
For both "cartesian" and "cylindrical" distributions, this element
specifies the distribution of z-coordinates. The necessary
sub-elements/attributes are those of a univariate probability
distribution (see the description in :ref:`univariate`).
:r:
For a "spherical" distribution, this element specifies the distribution
of r-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
For "cylindrical" and "spherical" distributions, this element specifies
the distribution of r-coordinates (cylindrical radius and spherical
radius, respectively). The necessary sub-elements/attributes are those
of a univariate probability distribution (see the description in
:ref:`univariate`).
:theta:
@ -498,14 +505,14 @@ attributes/sub-elements:
:ref:`univariate`).
:phi:
For a "spherical" distribution, this element specifies the distribution
of phi-coordinates. The necessary sub-elements/attributes are those of a
univariate probability distribution (see the description in
:ref:`univariate`).
For "cylindrical" and "spherical" distributions, this element specifies
the distribution of phi-coordinates. The necessary
sub-elements/attributes are those of a univariate probability
distribution (see the description in :ref:`univariate`).
:origin:
For a "spherical" distribution, this element specifies the coordinates of
the center of the sphere.
For "cylindrical and "spherical" distributions, this element specifies
the coordinates for the origin of the coordinate system.
:angle:
An element specifying the angular distribution of source sites. This element

View file

@ -4,7 +4,7 @@
License Agreement
=================
Copyright © 2011-2019 Massachusetts Institute of Technology and OpenMC contributors
Copyright © 2011-2020 Massachusetts Institute of Technology and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

View file

@ -46,6 +46,7 @@ Spatial Distributions
openmc.stats.Spatial
openmc.stats.CartesianIndependent
openmc.stats.CylindricalIndependent
openmc.stats.SphericalIndependent
openmc.stats.Box
openmc.stats.Point

View file

@ -140,7 +140,7 @@ extern "C" {
const int* indices, int n_elements,
int dim, double spectral,
const int* cmfd_indices,
const int* map);
const int* map, bool use_all_threads);
//! Runs a Gauss Seidel linear solver to solve CMFD matrix equations
//! linear solver

View file

@ -38,6 +38,26 @@ private:
UPtrDist z_; //!< Distribution of z coordinates
};
//==============================================================================
//! Distribution of points specified by cylindrical coordinates r,phi,z
//==============================================================================
class CylindricalIndependent : public SpatialDistribution {
public:
explicit CylindricalIndependent(pugi::xml_node node);
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const;
private:
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist phi_; //!< Distribution of phi coordinates
UPtrDist z_; //!< Distribution of z coordinates
Position origin_; //!< Cartesian coordinates of the cylinder center
};
//==============================================================================
//! Distribution of points specified by spherical coordinates r,theta,phi
//==============================================================================

View file

@ -54,7 +54,7 @@ Indicates the default path to an HDF5 file that contains multi-group cross
section libraries if the user has not specified the <cross_sections> tag in
.I materials.xml\fP.
.SH LICENSE
Copyright \(co 2011-2019 Massachusetts Institute of Technology and OpenMC
Copyright \(co 2011-2020 Massachusetts Institute of Technology and OpenMC
contributors.
.PP
Permission is hereby granted, free of charge, to any person obtaining a copy of

View file

@ -196,13 +196,10 @@ class CMFDRun(object):
----------
tally_begin : int
Batch number at which CMFD tallies should begin accummulating
feedback_begin: int
Batch number at which CMFD feedback should be turned on
solver_begin: int
Batch number at which CMFD solver should start executing
ref_d : list of floats
List of reference diffusion coefficients to fix CMFD parameters to
dhat_reset : bool
Indicate whether :math:`\widehat{D}` nonlinear CMFD parameters should
be reset to zero before solving CMFD eigenproblem.
display : dict
Dictionary indicating which CMFD results to output. Note that CMFD
k-effective will always be outputted. Acceptable keys are:
@ -298,6 +295,8 @@ class CMFDRun(object):
Time for building CMFD matrices, in seconds
time_cmfdsolve : float
Time for solving CMFD matrix equations, in seconds
use_all_threads : bool
Whether to use all threads allocated to OpenMC for CMFD solver
intracomm : mpi4py.MPI.Intracomm or None
MPI intercommunicator for running MPI commands
@ -310,9 +309,8 @@ class CMFDRun(object):
"""
# Variables that users can modify
self._tally_begin = 1
self._feedback_begin = 1
self._ref_d = []
self._dhat_reset = False
self._solver_begin = 1
self._ref_d = np.array([])
self._display = {'balance': False, 'dominance': False,
'entropy': False, 'source': False}
self._downscatter = False
@ -332,6 +330,7 @@ class CMFDRun(object):
self._window_type = 'none'
self._window_size = 10
self._intracomm = None
self._use_all_threads = False
# External variables used during runtime but users cannot control
self._set_reference_params = False
@ -339,7 +338,6 @@ class CMFDRun(object):
self._egrid = None
self._albedo = None
self._coremap = None
self._n_resets = 0
self._mesh_id = None
self._tally_ids = None
self._energy_filters = None
@ -418,17 +416,13 @@ class CMFDRun(object):
return self._tally_begin
@property
def feedback_begin(self):
return self._feedback_begin
def solver_begin(self):
return self._solver_begin
@property
def ref_d(self):
return self._ref_d
@property
def dhat_reset(self):
return self._dhat_reset
@property
def display(self):
return self._display
@ -501,6 +495,10 @@ class CMFDRun(object):
def indices(self):
return self._indices
@property
def use_all_threads(self):
return self._use_all_threads
@property
def cmfd_src(self):
return self._cmfd_src
@ -531,11 +529,12 @@ class CMFDRun(object):
check_greater_than('CMFD tally begin batch', begin, 0)
self._tally_begin = begin
@feedback_begin.setter
def feedback_begin(self, begin):
@solver_begin.setter
def solver_begin(self, begin):
check_type('CMFD feedback begin batch', begin, Integral)
check_greater_than('CMFD feedback begin batch', begin, 0)
self._feedback_begin = begin
self._solver_begin = begin
@ref_d.setter
def ref_d(self, diff_params):
@ -543,11 +542,6 @@ class CMFDRun(object):
Iterable, Real)
self._ref_d = np.array(diff_params)
@dhat_reset.setter
def dhat_reset(self, dhat_reset):
check_type('CMFD Dhat reset', dhat_reset, bool)
self._dhat_reset = dhat_reset
@display.setter
def display(self, display):
check_type('display', display, Mapping)
@ -691,6 +685,11 @@ class CMFDRun(object):
check_length('Gauss-Seidel tolerance', gauss_seidel_tolerance, 2)
self._gauss_seidel_tolerance = gauss_seidel_tolerance
@use_all_threads.setter
def use_all_threads(self, use_all_threads):
check_type('CMFD use all threads', use_all_threads, bool)
self._use_all_threads = use_all_threads
def run(self, **kwargs):
"""Run OpenMC with coarse mesh finite difference acceleration
@ -762,22 +761,23 @@ class CMFDRun(object):
calling :func:`openmc.lib.simulation_init`
"""
# Configure CMFD parameters and tallies
# Configure CMFD parameters
self._configure_cmfd()
# Initialize all arrays used for CMFD solver
self._allocate_cmfd()
# Create tally objects
self._create_cmfd_tally()
# Compute and store array indices used to build cross section
# arrays
self._precompute_array_indices()
if openmc.lib.master():
# Compute and store array indices used to build cross section
# arrays
self._precompute_array_indices()
# Compute and store row and column indices used to build CMFD
# matrices
self._precompute_matrix_indices()
# Compute and store row and column indices used to build CMFD
# matrices
self._precompute_matrix_indices()
# Initialize all variables used for linear solver in C++
self._initialize_linsolver()
# Initialize all variables used for linear solver in C++
self._initialize_linsolver()
# Initialize simulation
openmc.lib.simulation_init()
@ -801,13 +801,8 @@ class CMFDRun(object):
# Run next batch
status = openmc.lib.next_batch()
# Perform CMFD calculation if on
if self._cmfd_on:
self._execute_cmfd()
# Write CMFD output if CMFD on for current batch
if openmc.lib.master():
self._write_cmfd_output()
# Perform CMFD calculations
self._execute_cmfd()
# Write CMFD data to statepoint
if openmc.lib.is_statepoint_batch():
@ -823,8 +818,9 @@ class CMFDRun(object):
# Finalize simuation
openmc.lib.simulation_finalize()
# Print out CMFD timing statistics
self._write_cmfd_timing_stats()
if openmc.lib.master():
# Print out CMFD timing statistics
self._write_cmfd_timing_stats()
def statepoint_write(self, filename=None):
"""Write all simulation parameters to statepoint
@ -865,7 +861,7 @@ class CMFDRun(object):
cmfd_group = f.create_group("cmfd")
cmfd_group.attrs['cmfd_on'] = self._cmfd_on
cmfd_group.attrs['feedback'] = self._feedback
cmfd_group.attrs['feedback_begin'] = self._feedback_begin
cmfd_group.attrs['solver_begin'] = self._solver_begin
cmfd_group.attrs['mesh_id'] = self._mesh_id
cmfd_group.attrs['tally_begin'] = self._tally_begin
cmfd_group.attrs['time_cmfd'] = self._time_cmfd
@ -921,7 +917,7 @@ class CMFDRun(object):
args = temp_loss.indptr, len(temp_loss.indptr), \
temp_loss.indices, len(temp_loss.indices), n, \
self._spectral, self._indices, coremap
self._spectral, self._indices, coremap, self._use_all_threads
return openmc.lib._dll.openmc_initialize_linsolver(*args)
def _write_cmfd_output(self):
@ -948,53 +944,33 @@ class CMFDRun(object):
def _write_cmfd_timing_stats(self):
"""Write CMFD timing stats to buffer after finalizing simulation"""
if openmc.lib.master():
outstr = ("=====================> "
"CMFD TIMING STATISTICS <====================\n\n"
" Time in CMFD = {:.5E} seconds\n"
" Building matrices = {:.5E} seconds\n"
" Solving matrices = {:.5E} seconds\n")
print(outstr.format(self._time_cmfd, self._time_cmfdbuild,
self._time_cmfdsolve))
sys.stdout.flush()
outstr = ("=====================> "
"CMFD TIMING STATISTICS <====================\n\n"
" Time in CMFD = {:.5e} seconds\n"
" Building matrices = {:.5e} seconds\n"
" Solving matrices = {:.5e} seconds\n")
print(outstr.format(self._time_cmfd, self._time_cmfdbuild,
self._time_cmfdsolve))
sys.stdout.flush()
def _configure_cmfd(self):
"""Initialize CMFD parameters and set CMFD input variables"""
# Check if restarting simulation from statepoint file
if not openmc.lib.settings.restart_run:
# Read in cmfd input defined in Python
self._read_cmfd_input()
# Set up CMFD coremap
self._set_coremap()
# Extract spatial and energy indices
nx, ny, nz, ng = self._indices
# Allocate parameters that need to stored for tally window
self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0))
self._flux_rate = np.zeros((nx, ny, nz, ng, 0))
self._total_rate = np.zeros((nx, ny, nz, ng, 0))
self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0))
self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0))
self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0))
self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0))
# Initialize timers
self._time_cmfd = 0.0
self._time_cmfdbuild = 0.0
self._time_cmfdsolve = 0.0
# Initialize parameters for CMFD tally windows
self._set_tally_window()
# Define all variables necessary for running CMFD
self._initialize_cmfd()
else:
# Reset CMFD parameters from statepoint file
path_statepoint = openmc.lib.settings.path_statepoint
self._reset_cmfd(path_statepoint)
def _read_cmfd_input(self):
"""Sets values of additional instance variables based on user input"""
def _initialize_cmfd(self):
"""Sets values of CMFD instance variables based on user input,
separating between variables that only exist on all processes
and those that only exist on the master process
"""
# Print message to user and flush output to stdout
if openmc.lib.settings.verbosity >= 7 and openmc.lib.master():
print(' Configuring CMFD parameters for simulation')
@ -1024,34 +1000,55 @@ class CMFDRun(object):
self._indices[3] = 1
self._energy_filters = False
# Set global albedo
if self._mesh.albedo is not None:
self._albedo = np.array(self._mesh.albedo)
else:
self._albedo = np.array([1., 1., 1., 1., 1., 1.])
# Get acceleration map, otherwise set all regions to be accelerated
if self._mesh.map is not None:
check_length('CMFD coremap', self._mesh.map,
np.product(self._indices[0:3]))
self._coremap = np.array(self._mesh.map)
if openmc.lib.master():
self._coremap = np.array(self._mesh.map)
else:
self._coremap = np.ones((np.product(self._indices[0:3])),
dtype=int)
if openmc.lib.master():
self._coremap = np.ones((np.product(self._indices[0:3])),
dtype=int)
# Check CMFD tallies accummulated before feedback turned on
if self._feedback and self._feedback_begin < self._tally_begin:
if self._feedback and self._solver_begin < self._tally_begin:
raise ValueError('Tally begin must be less than or equal to '
'feedback begin')
'CMFD begin')
# Set number of batches where cmfd tallies should be reset
self._n_resets = len(self._reset)
# Initialize parameters for CMFD tally windows
self._set_tally_window()
# Create tally objects
self._create_cmfd_tally()
# Define all variables that will exist only on master process
if openmc.lib.master():
# Set global albedo
if self._mesh.albedo is not None:
self._albedo = np.array(self._mesh.albedo)
else:
self._albedo = np.array([1., 1., 1., 1., 1., 1.])
# Set up CMFD coremap
self._set_coremap()
# Extract spatial and energy indices
nx, ny, nz, ng = self._indices
# Allocate parameters that need to be stored for tally window
self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0))
self._flux_rate = np.zeros((nx, ny, nz, ng, 0))
self._total_rate = np.zeros((nx, ny, nz, ng, 0))
self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0))
self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0))
self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0))
self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0))
# Initialize timers
self._time_cmfd = 0.0
self._time_cmfdbuild = 0.0
self._time_cmfdsolve = 0.0
def _reset_cmfd(self, filename):
"""Reset all CMFD parameters from statepoint
"""Reset all CMFD parameters from statepoint
Parameters
----------
@ -1070,32 +1067,28 @@ class CMFDRun(object):
print(' Loading CMFD data from {}...'.format(filename))
sys.stdout.flush()
cmfd_group = f['cmfd']
# Define variables that exist on all processes
self._cmfd_on = cmfd_group.attrs['cmfd_on']
self._feedback = cmfd_group.attrs['feedback']
self._feedback_begin = cmfd_group.attrs['feedback_begin']
self._solver_begin = cmfd_group.attrs['solver_begin']
self._tally_begin = cmfd_group.attrs['tally_begin']
self._time_cmfd = cmfd_group.attrs['time_cmfd']
self._time_cmfdbuild = cmfd_group.attrs['time_cmfdbuild']
self._time_cmfdsolve = cmfd_group.attrs['time_cmfdsolve']
self._window_size = cmfd_group.attrs['window_size']
self._window_type = cmfd_group.attrs['window_type']
self._k_cmfd = list(cmfd_group['k_cmfd'])
self._dom = list(cmfd_group['dom'])
self._src_cmp = list(cmfd_group['src_cmp'])
self._balance = list(cmfd_group['balance'])
self._entropy = list(cmfd_group['entropy'])
self._reset = list(cmfd_group['reset'])
self._albedo = cmfd_group['albedo'][()]
self._coremap = cmfd_group['coremap'][()]
self._egrid = cmfd_group['egrid'][()]
self._indices = cmfd_group['indices'][()]
self._current_rate = cmfd_group['current_rate'][()]
self._flux_rate = cmfd_group['flux_rate'][()]
self._nfiss_rate = cmfd_group['nfiss_rate'][()]
self._openmc_src_rate = cmfd_group['openmc_src_rate'][()]
self._p1scatt_rate = cmfd_group['p1scatt_rate'][()]
self._scatt_rate = cmfd_group['scatt_rate'][()]
self._total_rate = cmfd_group['total_rate'][()]
default_egrid = np.array([_ENERGY_MIN_NEUTRON,
_ENERGY_MAX_NEUTRON])
self._energy_filters = not np.array_equal(self._egrid,
default_egrid)
self._window_size = cmfd_group.attrs['window_size']
self._window_type = cmfd_group.attrs['window_type']
self._reset_every = (self._window_type == 'expanding' or
self._window_type == 'rolling')
# Overwrite CMFD mesh properties
cmfd_mesh_name = 'mesh ' + str(cmfd_group.attrs['mesh_id'])
@ -1105,49 +1098,21 @@ class CMFDRun(object):
self._mesh.upper_right = cmfd_mesh['upper_right'][()]
self._mesh.width = cmfd_mesh['width'][()]
# Store tally ids from statepoint run
sp_tally_ids = list(cmfd_group['tally_ids'])
# Set CMFD variables not in statepoint file
default_egrid = np.array([_ENERGY_MIN_NEUTRON, _ENERGY_MAX_NEUTRON])
self._energy_filters = not np.array_equal(self._egrid, default_egrid)
self._n_resets = len(self._reset)
self._mat_dim = np.max(self._coremap) + 1
self._reset_every = (self._window_type == 'expanding' or
self._window_type == 'rolling')
# Recreate CMFD tallies in memory
self._create_cmfd_tally()
def _allocate_cmfd(self):
"""Allocates all numpy arrays and lists used in CMFD algorithm"""
# Extract spatial and energy indices
nx, ny, nz, ng = self._indices
# Allocate dimensions for each mesh cell
self._hxyz = np.zeros((nx, ny, nz, 3))
self._hxyz[:] = openmc.lib.meshes[self._mesh_id].width
# Allocate flux, cross sections and diffusion coefficient
self._flux = np.zeros((nx, ny, nz, ng))
self._totalxs = np.zeros((nx, ny, nz, ng))
self._p1scattxs = np.zeros((nx, ny, nz, ng))
self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing
self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing
self._diffcof = np.zeros((nx, ny, nz, ng))
# Allocate dtilde and dhat
self._dtilde = np.zeros((nx, ny, nz, ng, 6))
self._dhat = np.zeros((nx, ny, nz, ng, 6))
# Set reference diffusion parameters
if self._ref_d:
self._set_reference_params = True
# Check length of reference diffusion parameters equal to number of
# energy groups
if len(self._ref_d) != self._indices[3]:
raise OpenMCError('Number of reference diffusion parameters '
'must equal number of CMFD energy groups')
# Define variables that exist only on master process
if openmc.lib.master():
self._time_cmfd = cmfd_group.attrs['time_cmfd']
self._time_cmfdbuild = cmfd_group.attrs['time_cmfdbuild']
self._time_cmfdsolve = cmfd_group.attrs['time_cmfdsolve']
self._albedo = cmfd_group['albedo'][()]
self._coremap = cmfd_group['coremap'][()]
self._current_rate = cmfd_group['current_rate'][()]
self._flux_rate = cmfd_group['flux_rate'][()]
self._nfiss_rate = cmfd_group['nfiss_rate'][()]
self._openmc_src_rate = cmfd_group['openmc_src_rate'][()]
self._p1scatt_rate = cmfd_group['p1scatt_rate'][()]
self._scatt_rate = cmfd_group['scatt_rate'][()]
self._total_rate = cmfd_group['total_rate'][()]
self._mat_dim = np.max(self._coremap) + 1
def _set_tally_window(self):
"""Sets parameters to handle different tally window options"""
@ -1169,47 +1134,57 @@ class CMFDRun(object):
# Add 1 as next_batch has not been called yet
current_batch = openmc.lib.current_batch() + 1
# Check to activate CMFD diffusion and possible feedback
# Check to activate CMFD tallies
if self._tally_begin == current_batch:
# Check to activate CMFD solver and possible feedback
if self._solver_begin == current_batch:
self._cmfd_on = True
# Check to reset tallies
if ((self._n_resets > 0 and current_batch in self._reset)
if ((len(self._reset) > 0 and current_batch in self._reset)
or self._reset_every):
self._cmfd_tally_reset()
def _execute_cmfd(self):
"""Runs CMFD calculation on master node"""
# Run CMFD on single processor on master
if openmc.lib.master():
# Start CMFD timer
time_start_cmfd = time.time()
# Create CMFD data from OpenMC tallies
self._set_up_cmfd()
if openmc.lib.current_batch() >= self._tally_begin:
# Calculate all cross sections based on tally window averages
self._compute_xs()
# Call solver
self._cmfd_solver_execute()
# Execute CMFD algorithm if CMFD on for current batch
if self._cmfd_on:
# Run CMFD on single processor on master
if openmc.lib.master():
# Create CMFD data based on OpenMC tallies
self._set_up_cmfd()
# Store k-effective
self._k_cmfd.append(self._keff)
# Call solver
self._cmfd_solver_execute()
# Check to perform adjoint on last batch
if (openmc.lib.current_batch() == openmc.lib.settings.batches
and self._run_adjoint):
self._cmfd_solver_execute(adjoint=True)
# Store k-effective
self._k_cmfd.append(self._keff)
# Calculate fission source
self._calc_fission_source()
# Check to perform adjoint on last batch
if (openmc.lib.current_batch() == openmc.lib.settings.batches
and self._run_adjoint):
self._cmfd_solver_execute(adjoint=True)
# Calculate weight factors
self._cmfd_reweight(True)
# Calculate fission source
self._calc_fission_source()
# Calculate weight factors
self._cmfd_reweight()
# Stop CMFD timer
if openmc.lib.master():
time_stop_cmfd = time.time()
self._time_cmfd += time_stop_cmfd - time_start_cmfd
if self._cmfd_on:
# Write CMFD output if CMFD on for current batch
self._write_cmfd_output()
def _cmfd_tally_reset(self):
"""Resets all CMFD tallies in memory"""
@ -1228,9 +1203,6 @@ class CMFDRun(object):
"""Configures CMFD object for a CMFD eigenvalue calculation
"""
# Calculate all cross sections based on tally window averages
self._compute_xs()
# Compute effective downscatter cross section
if self._downscatter:
self._compute_effective_downscatter()
@ -1422,96 +1394,85 @@ class CMFDRun(object):
self._src_cmp.append(np.sqrt(1.0 / self._norm
* np.sum((self._cmfd_src - self._openmc_src)**2)))
def _cmfd_reweight(self, new_weights):
"""Performs weighting of particles in source bank
def _cmfd_reweight(self):
"""Performs weighting of particles in source bank"""
# Get spatial dimensions and energy groups
nx, ny, nz, ng = self._indices
Parameters
----------
new_weights : bool
Whether to reweight particles or not
# Count bank site in mesh and reverse due to egrid structured
outside = self._count_bank_sites()
"""
# Compute new weight factors
if new_weights:
# 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')
# Get spatial dimensions and energy groups
nx, ny, nz, ng = self._indices
# 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)
# Count bank site in mesh and reverse due to egrid structured
outside = self._count_bank_sites()
# Define target reshape dimensions for sourcecounts. This
# defines how self._sourcecounts is ordered by dimension
target_shape = [nz, ny, nx, ng]
# 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')
# 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)
# 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)
# Flip index of energy dimension
sourcecounts = np.flip(sourcecounts, axis=3)
# Define target reshape dimensions for sourcecounts. This
# defines how self._sourcecounts is ordered by dimension
target_shape = [nz, ny, nx, ng]
# 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))
# 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)
if not self._feedback:
return
# Flip index of energy dimension
sourcecounts = np.flip(sourcecounts, axis=3)
# Broadcast weight factors to all procs
if have_mpi:
self._weightfactors = self._intracomm.bcast(
self._weightfactors)
# 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))
m = openmc.lib.meshes[self._mesh_id]
energy = self._egrid
ng = self._indices[3]
if (not self._feedback
or openmc.lib.current_batch() < self._feedback_begin):
return
# Get locations and energies of all particles in source bank
source_xyz = openmc.lib.source_bank()['r']
source_energies = openmc.lib.source_bank()['E']
# Broadcast weight factors to all procs
if have_mpi:
self._weightfactors = self._intracomm.bcast(
self._weightfactors)
# Convert xyz location to the CMFD mesh index
mesh_ijk = np.floor((source_xyz - m.lower_left)/m.width).astype(int)
m = openmc.lib.meshes[self._mesh_id]
energy = self._egrid
ng = self._indices[3]
# 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)
# Get locations and energies of all particles in source bank
source_xyz = openmc.lib.source_bank()['r']
source_energies = openmc.lib.source_bank()['E']
# 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]
# 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, 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 pt below energy grid')
sys.stdout.flush()
if openmc.lib.master() and np.any(source_energies > energy[-1]):
print(' WARNING: Source pt above energy grid')
sys.stdout.flush()
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
@ -1556,7 +1517,7 @@ class CMFDRun(object):
energy_bins[idx] = ng - 1
idx = np.where((source_energies >= energy[0]) &
(source_energies <= energy[-1]))
energy_bins[idx] = np.digitize(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
@ -1866,8 +1827,8 @@ class CMFDRun(object):
if self._power_monitor and openmc.lib.master():
str1 = ' {:d}:'.format(iter)
str2 = 'k-eff: {:0.8f}'.format(k_n)
str3 = 'k-error: {:.5E}'.format(kerr)
str4 = 'src-error: {:.5E}'.format(serr)
str3 = 'k-error: {:.5e}'.format(kerr)
str4 = 'src-error: {:.5e}'.format(serr)
str5 = ' {:d}'.format(innerits)
print('{:8s}{:20s}{:25s}{:s}{:s}'.format(str1, str2, str3, str4,
str5))
@ -1927,33 +1888,12 @@ class CMFDRun(object):
# Get tallies in-memory
tallies = openmc.lib.tallies
# Ravel coremap as 1d array similar to how tally data is arranged
coremap = np.ravel(self._coremap.swapaxes(0, 2))
# Set conditional numpy array as boolean vector based on coremap
# Repeat each value for number of groups in problem
is_cmfd_accel = np.repeat(coremap != _CMFD_NOACCEL, ng)
is_accel = self._coremap != _CMFD_NOACCEL
# Get flux from CMFD tally 0
tally_id = self._tally_ids[0]
tally_results = tallies[tally_id].results[:,0,1]
flux = np.where(is_cmfd_accel, tally_results, 0.)
# Detect zero flux, abort if located
if np.any(flux[is_cmfd_accel] < _TINY_BIT):
# Get index of zero flux in flux array
idx = np.argmax(np.where(is_cmfd_accel, flux, 1) < _TINY_BIT)
# Convert scalar idx to index in flux matrix
mat_idx = np.unravel_index(idx, self._flux.shape)
# Throw error message (one-based indexing)
# Index of group is flipped
err_message = 'Detected zero flux without coremap overlay' + \
' at mesh: (' + \
', '.join(str(i+1) for i in mat_idx[:-1]) + \
') in group ' + str(ng-mat_idx[-1])
raise OpenMCError(err_message)
flux = tallies[tally_id].results[:,0,1]
# Define target tally reshape dimensions. This defines how openmc
# tallies are ordered by dimension
@ -1971,7 +1911,23 @@ class CMFDRun(object):
self._flux_rate = np.append(self._flux_rate, reshape_flux, axis=4)
# Compute flux as aggregate of banked flux_rate over tally window
self._flux = np.sum(self._flux_rate, axis=4)
self._flux = np.where(is_accel[..., np.newaxis],
np.sum(self._flux_rate, axis=4), 0.0)
# Detect zero flux, abort if located and cmfd is on
zero_flux = np.logical_and(self._flux < _TINY_BIT,
is_accel[..., np.newaxis])
if np.any(zero_flux) and self._cmfd_on:
# Get index of first zero flux in flux array
idx = np.argwhere(zero_flux)[0]
# Throw error message (one-based indexing)
# Index of group is flipped
err_message = 'Detected zero flux without coremap overlay' + \
' at mesh: (' + \
', '.join(str(i+1) for i in idx[:-1]) + \
') in group ' + str(ng-idx[-1])
raise OpenMCError(err_message)
# Get total rr from CMFD tally 0
totalrr = tallies[tally_id].results[:,1,1]
@ -2072,10 +2028,7 @@ class CMFDRun(object):
# Get surface currents from CMFD tally 2
tally_id = self._tally_ids[2]
tally_results = tallies[tally_id].results[:,0,1]
# Filter tally results to include only accelerated regions
current = np.where(np.repeat(is_cmfd_accel, 12), tally_results, 0.)
current = tallies[tally_id].results[:,0,1]
# Define target tally reshape dimensions for current
target_tally_shape = [nz, ny, nx, 12, ng, 1]
@ -2094,7 +2047,8 @@ class CMFDRun(object):
axis=5)
# Compute current as aggregate of banked current_rate over tally window
self._current = np.sum(self._current_rate, axis=5)
self._current = np.where(is_accel[..., np.newaxis, np.newaxis],
np.sum(self._current_rate, axis=5), 0.0)
# Get p1 scatter rr from CMFD tally 3
tally_id = self._tally_ids[3]
@ -2202,19 +2156,20 @@ class CMFDRun(object):
# Compute scattering rr by broadcasting flux in outgoing energy and
# summing over incoming energy
scattering = np.sum(self._scattxs * self._flux[:,:,:,:,np.newaxis],
scattering = np.sum(self._scattxs * self._flux[:,:,:,:, np.newaxis],
axis=3)
# Compute fission rr by broadcasting flux in outgoing energy and
# summing over incoming energy
fission = np.sum(self._nfissxs * self._flux[:,:,:,:,np.newaxis],
fission = np.sum(self._nfissxs * self._flux[:,:,:,:, np.newaxis],
axis=3)
# Compute residual
res = leakage + interactions - scattering - (1.0 / keff) * fission
# Normalize res by flux and bank res
self._resnb = np.divide(res, self._flux, where=self._flux > 0)
self._resnb = np.divide(res, self._flux, where=self._flux > 0,
out=np.zeros_like(self._flux))
# Calculate RMS and record for this batch
self._balance.append(np.sqrt(
@ -2222,12 +2177,37 @@ class CMFDRun(object):
(ng * num_accel)))
def _precompute_array_indices(self):
"""Computes the indices used to populate certain cross section arrays.
These indices are used in _compute_dtilde and _compute_dhat
"""Initializes cross section arrays and computes the indices
used to populate dtilde and dhat
"""
# Extract spatial indices
nx, ny, nz = self._indices[:3]
nx, ny, nz, ng = self._indices
# Allocate dimensions for each mesh cell
self._hxyz = np.zeros((nx, ny, nz, 3))
self._hxyz[:] = openmc.lib.meshes[self._mesh_id].width
# Allocate flux, cross sections and diffusion coefficient
self._flux = np.zeros((nx, ny, nz, ng))
self._totalxs = np.zeros((nx, ny, nz, ng))
self._p1scattxs = np.zeros((nx, ny, nz, ng))
self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing
self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing
self._diffcof = np.zeros((nx, ny, nz, ng))
# Allocate dtilde and dhat
self._dtilde = np.zeros((nx, ny, nz, ng, 6))
self._dhat = np.zeros((nx, ny, nz, ng, 6))
# Set reference diffusion parameters
if self._ref_d.size > 0:
self._set_reference_params = True
# Check length of reference diffusion parameters equal to number of
# energy groups
if self._ref_d.size != self._indices[3]:
raise OpenMCError('Number of reference diffusion parameters '
'must equal number of CMFD energy groups')
# Logical for determining whether region of interest is accelerated
# region

View file

@ -23,13 +23,16 @@ class DataLibrary(EqualityMixin):
def __init__(self):
self.libraries = []
def get_by_material(self, name):
def get_by_material(self, name, data_type='neutron'):
"""Return the library dictionary containing a given material.
Parameters
----------
name : str
Name of material, e.g. 'Am241'
data_type : str
Name of data type, e.g. 'neutron', 'photon', 'wmp',
or 'thermal'
Returns
-------
@ -39,7 +42,7 @@ class DataLibrary(EqualityMixin):
"""
for library in self.libraries:
if name in library['materials']:
if name in library['materials'] and data_type in library['type']:
return library
return None

View file

@ -47,7 +47,7 @@ _dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
_dll.openmc_get_keff.restype = c_int
_dll.openmc_get_keff.errcheck = _error_handler
_init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int,
c_double, _array_1d_int, _array_1d_int]
c_double, _array_1d_int, _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

View file

@ -999,6 +999,11 @@ class Library(object):
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
elif 'transport' in self.mgxs_types and self.correction == 'P0':
mymgxs = self.get_mgxs(domain, 'transport')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
subdomain=subdomain)
elif 'total' in self.mgxs_types:
mymgxs = self.get_mgxs(domain, 'total')
xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
@ -1167,6 +1172,16 @@ class Library(object):
np.subtract(xsdata._total[i], np.sum(np.sum(
xsdata._scatter_matrix[i][:, :, :, :, :],
axis=4), axis=3))
# if only scatter matrices have been tallied, multiplicity cannot
# be accounted for
else:
msg = 'Scatter multiplicity (such as (n,xn) reactions) '\
'are ignored since multiplicity or nu-scatter matrices '\
'were not tallied for ' + xsdata_name
warn(msg, RuntimeWarning)
xsdata.set_scatter_matrix_mgxs(scatt_mgxs, xs_type=xs_type,
nuclide=[nuclide],
subdomain=subdomain)
return xsdata
@ -1410,6 +1425,9 @@ class Library(object):
fixed source problem could be the target.
- Fission and kappa-fission are not required as they are only
needed to support tallies the user may wish to request.
- Scattering multiplicity should have been tallied for increased model
accuracy, either using a multiplicity or scatter and nu-scatter matrix
tally.
See also
--------
@ -1460,6 +1478,15 @@ class Library(object):
'"scatter matrix", or "consistent scatter matrix" MGXS '
'type is required.')
# Make sure there is some kind of a scattering multiplicity matrix data
if 'multiplicity matrix' not in self.mgxs_types and \
('scatter matrix' not in self.mgxs_types or
'nu-scatter matrix' not in self.mgxs_types) and\
('consistent scatter matrix' not in self.mgxs_types or
'consistent nu-scatter matrix' not in self.mgxs_types):
warn('A "multiplicity matrix" or both a "scatter" and "nu-scatter" '
'matrix MGXS type(s) should be provided.')
# Ensure absorption is present
if 'absorption' not in self.mgxs_types:
error_flag = True

View file

@ -1638,9 +1638,9 @@ class XSdata(object):
xs_type=xs_type, moment=0,
subdomains=subdomain)
if scatter.scatter_format == 'histogram':
scatt = np.sum(scatt, axis=0)
scatt = np.sum(scatt, axis=2)
if nuscatter.scatter_format == 'histogram':
nuscatt = np.sum(nuscatt, axis=0)
nuscatt = np.sum(nuscatt, axis=2)
self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
self._multiplicity_matrix[i] = \

View file

@ -562,14 +562,16 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
if isinstance(this, openmc.Material):
for sab_name in this._sab:
sab = openmc.data.ThermalScattering.from_hdf5(
library.get_by_material(sab_name)['path'])
library.get_by_material(sab_name, data_type='thermal')['path'])
for nuc in sab.nuclides:
sabs[nuc] = library.get_by_material(sab_name)['path']
sabs[nuc] = library.get_by_material(sab_name,
data_type='thermal')['path']
else:
if sab_name:
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
for nuc in sab.nuclides:
sabs[nuc] = library.get_by_material(sab_name)['path']
sabs[nuc] = library.get_by_material(sab_name,
data_type='thermal')['path']
# Now we can create the data sets to be plotted
xs = {}

View file

@ -272,6 +272,8 @@ class Spatial(metaclass=ABCMeta):
distribution = get_text(elem, 'type')
if distribution == 'cartesian':
return CartesianIndependent.from_xml_element(elem)
elif distribution == 'cylindrical':
return CylindricalIndependent.from_xml_element(elem)
elif distribution == 'spherical':
return SphericalIndependent.from_xml_element(elem)
elif distribution == 'box' or distribution == 'fission':
@ -377,7 +379,7 @@ class CartesianIndependent(Spatial):
class SphericalIndependent(Spatial):
"""Spatial distribution represented in spherical coordinates.
r"""Spatial distribution represented in spherical coordinates.
This distribution allows one to specify coordinates whose :math:`r`,
:math:`\theta`, and :math:`\phi` components are sampled independently from
@ -386,26 +388,31 @@ class SphericalIndependent(Spatial):
Parameters
----------
r : openmc.stats.Univariate
Distribution of r-coordinates
Distribution of r-coordinates in a reference frame specified by
the origin parameter
theta : openmc.stats.Univariate
Distribution of theta-coordinates (angle relative to the z-axis)
Distribution of theta-coordinates (angle relative to the z-axis) in a
reference frame specified by the origin parameter
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle)
Distribution of phi-coordinates (azimuthal angle) in a reference frame
specified by the origin parameter
origin: Iterable of float, optional
coordinates (x0, y0, z0) of the center of the sphere. Defaults to
(0.0, 0.0, 0.0)
coordinates (x0, y0, z0) of the center of the spherical reference frame
for the source. Defaults to (0.0, 0.0, 0.0)
Attributes
----------
r : openmc.stats.Univariate
Distribution of r-coordinates
Distribution of r-coordinates in the local reference frame
theta : openmc.stats.Univariate
Distribution of theta-coordinates (angle relative to the z-axis)
Distribution of theta-coordinates (angle relative to the z-axis) in the
local reference frame
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle)
Distribution of phi-coordinates (azimuthal angle) in the local
reference frame
origin: Iterable of float, optional
coordinates (x0, y0, z0) of the center of the sphere. Defaults to
(0.0, 0.0, 0.0)
coordinates (x0, y0, z0) of the center of the spherical reference
frame. Defaults to (0.0, 0.0, 0.0)
"""
@ -491,6 +498,126 @@ class SphericalIndependent(Spatial):
origin = [float(x) for x in elem.get('origin').split()]
return cls(r, theta, phi, origin=origin)
class CylindricalIndependent(Spatial):
"""Spatial distribution represented in cylindrical coordinates.
This distribution allows one to specify coordinates whose :math:`r`,
:math:`\phi`, and :math:`z` components are sampled independently from
one another and in a reference frame whose origin is specified by the
coordinates (x0, y0, z0).
Parameters
----------
r : openmc.stats.Univariate
Distribution of r-coordinates in a reference frame specified by the
origin parameter
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle) in a reference frame
specified by the origin parameter
z : openmc.stats.Univariate
Distribution of z-coordinates in a reference frame specified by the
origin parameter
origin: Iterable of float, optional
coordinates (x0, y0, z0) of the center of the cylindrical reference
frame. Defaults to (0.0, 0.0, 0.0)
Attributes
----------
r : openmc.stats.Univariate
Distribution of r-coordinates in the local reference frame
phi : openmc.stats.Univariate
Distribution of phi-coordinates (azimuthal angle) in the local
reference frame
z : openmc.stats.Univariate
Distribution of z-coordinates in the local reference frame
origin: Iterable of float, optional
coordinates (x0, y0, z0) of the center of the cylindrical reference
frame. Defaults to (0.0, 0.0, 0.0)
"""
def __init__(self, r, phi, z, origin=(0.0, 0.0, 0.0)):
super().__init__()
self.r = r
self.phi = phi
self.z = z
self.origin = origin
@property
def r(self):
return self._r
@property
def phi(self):
return self._phi
@property
def z(self):
return self._z
@property
def origin(self):
return self._origin
@r.setter
def r(self, r):
cv.check_type('r coordinate', r, Univariate)
self._r = r
@phi.setter
def phi(self, phi):
cv.check_type('phi coordinate', phi, Univariate)
self._phi = phi
@z.setter
def z(self, z):
cv.check_type('z coordinate', z, Univariate)
self._z = z
@origin.setter
def origin(self, origin):
cv.check_type('origin coordinates', origin, Iterable, Real)
origin = np.asarray(origin)
self._origin = origin
def to_xml_element(self):
"""Return XML representation of the spatial distribution
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing spatial distribution data
"""
element = ET.Element('space')
element.set('type', 'cylindrical')
element.append(self.r.to_xml_element('r'))
element.append(self.phi.to_xml_element('phi'))
element.append(self.z.to_xml_element('z'))
element.set("origin", ' '.join(map(str, self.origin)))
return element
@classmethod
def from_xml_element(cls, elem):
"""Generate spatial distribution from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.stats.CylindricalIndependent
Spatial distribution generated from XML element
"""
r = Univariate.from_xml_element(elem.find('r'))
phi = Univariate.from_xml_element(elem.find('phi'))
z = Univariate.from_xml_element(elem.find('z'))
origin = [float(x) for x in elem.get('origin').split()]
return cls(r, phi, z, origin=origin)
class Box(Spatial):
"""Uniform distribution of coordinates in a rectangular cuboid.

View file

@ -158,7 +158,7 @@ with tempfile.TemporaryDirectory() as tmpdir:
print('Creating compressed archive...')
test_tar = pwd / 'nndc_hdf5_test.tar.xz'
with tarfile.open(str(test_tar), 'w:xz') as txz:
txz.add('nndc_hdf5')
txz.add(output_dir)
# Change back to original directory
os.chdir(str(pwd))

View file

@ -3,6 +3,9 @@
#include <vector>
#include <cmath>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "xtensor/xtensor.hpp"
#include "openmc/error.h"
@ -29,6 +32,8 @@ int nx, ny, nz, ng;
xt::xtensor<int, 2> indexmap;
int use_all_threads;
} // namespace cmfd
//==============================================================================
@ -101,6 +106,7 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x,
for (int irb = 0; irb < 2; irb++) {
// Loop around matrix rows
#pragma omp parallel for reduction (+:err) if(cmfd::use_all_threads)
for (int irow = 0; irow < cmfd::dim; irow++) {
int g, i, j, k;
matrix_to_indices(irow, g, i, j, k);
@ -167,6 +173,7 @@ int cmfd_linsolver_2g(const double* A_data, const double* b, double* x,
for (int irb = 0; irb < 2; irb++) {
// Loop around matrix rows
#pragma omp parallel for reduction (+:err) if(cmfd::use_all_threads)
for (int irow = 0; irow < cmfd::dim; irow+=2) {
int g, i, j, k;
matrix_to_indices(irow, g, i, j, k);
@ -302,7 +309,7 @@ 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)
const int* map, bool use_all_threads)
{
// Store elements of indptr
for (int i = 0; i < len_indptr; i++)
@ -329,6 +336,9 @@ void openmc_initialize_linsolver(const int* indptr, int len_indptr,
cmfd::indexmap.resize({static_cast<size_t>(dim), 3});
set_indexmap(map);
}
// Use all threads allocated to OpenMC simulation to run CMFD solver
cmfd::use_all_threads = use_all_threads;
}
//==============================================================================

View file

@ -51,6 +51,72 @@ Position CartesianIndependent::sample(uint64_t* seed) const
return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};
}
//==============================================================================
// CylindricalIndependent implementation
//==============================================================================
CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
{
// Read distribution for r-coordinate
if (check_for_node(node, "r")) {
pugi::xml_node node_dist = node.child("r");
r_ = distribution_from_xml(node_dist);
} else {
// If no distribution was specified, default to a single point at r=0
double x[] {0.0};
double p[] {1.0};
r_ = std::make_unique<Discrete>(x, p, 1);
}
// Read distribution for phi-coordinate
if (check_for_node(node, "phi")) {
pugi::xml_node node_dist = node.child("phi");
phi_ = distribution_from_xml(node_dist);
} else {
// If no distribution was specified, default to a single point at phi=0
double x[] {0.0};
double p[] {1.0};
phi_ = std::make_unique<Discrete>(x, p, 1);
}
// Read distribution for z-coordinate
if (check_for_node(node, "z")) {
pugi::xml_node node_dist = node.child("z");
z_ = distribution_from_xml(node_dist);
} else {
// If no distribution was specified, default to a single point at z=0
double x[] {0.0};
double p[] {1.0};
z_ = std::make_unique<Discrete>(x, p, 1);
}
// Read cylinder center coordinates
if (check_for_node(node, "origin")) {
auto origin = get_node_array<double>(node, "origin");
if (origin.size() == 3) {
origin_ = origin;
} else {
std::stringstream err_msg;
err_msg << "Origin for cylindrical source distribution must be length 3";
fatal_error(err_msg);
}
} else {
// If no coordinates were specified, default to (0, 0, 0)
origin_ = {0.0, 0.0, 0.0};
}
}
Position CylindricalIndependent::sample(uint64_t* seed) const
{
double r = r_->sample(seed);
double phi = phi_->sample(seed);
double x = r*cos(phi) + origin_.x;
double y = r*sin(phi) + origin_.y;
double z = z_->sample(seed) + origin_.z;
return {x, y, z};
}
//==============================================================================
// SphericalIndependent implementation
//==============================================================================

View file

@ -44,7 +44,9 @@ void free_memory()
free_memory_mesh();
free_memory_tally();
free_memory_bank();
free_memory_cmfd();
if (mpi::master) {
free_memory_cmfd();
}
#ifdef DAGMC
free_memory_dagmc();
#endif

View file

@ -72,7 +72,7 @@ void title()
// Write version information
std::cout <<
" | The OpenMC Monte Carlo Code\n" <<
" Copyright | 2011-2019 MIT and OpenMC contributors\n" <<
" Copyright | 2011-2020 MIT and OpenMC contributors\n" <<
" License | http://openmc.readthedocs.io/en/latest/license.html\n" <<
" Version | " << VERSION_MAJOR << '.' << VERSION_MINOR << '.'
<< VERSION_RELEASE << (VERSION_DEV ? "-dev" : "") << '\n';

View file

@ -86,6 +86,8 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
type = get_node_value(node_space, "type", true, true);
if (type == "cartesian") {
space_ = UPtrSpace{new CartesianIndependent(node_space)};
} else if (type == "cylindrical") {
space_ = UPtrSpace{new CylindricalIndependent(node_space)};
} else if (type == "spherical") {
space_ = UPtrSpace{new SphericalIndependent(node_space)};
} else if (type == "box") {

View file

@ -24,7 +24,7 @@ def test_cmfd_physical_adjoint():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run_adjoint = True
@ -54,7 +54,7 @@ def test_cmfd_math_adjoint():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run_adjoint = True
@ -83,7 +83,7 @@ def test_cmfd_write_matrices():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
@ -131,7 +131,7 @@ def test_cmfd_feed():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
@ -140,3 +140,27 @@ def test_cmfd_feed():
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()
def test_cmfd_multithread():
"""Test 1 group CMFD solver with all available threads"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = (-10.0, -1.0, -1.0)
cmfd_mesh.upper_right = (10.0, 1.0, 1.0)
cmfd_mesh.dimension = (10, 1, 1)
cmfd_mesh.albedo = (0.0, 0.0, 1.0, 1.0, 1.0, 1.0)
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.use_all_threads = True
cmfd_run.run()
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()

View file

@ -17,7 +17,7 @@ def test_cmfd_feed_2g():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.downscatter = True

View file

@ -391,11 +391,6 @@ cmfd indices
1.000000E+00
1.000000E+00
k cmfd
1.143597E+00
1.163387E+00
1.173384E+00
1.171035E+00
1.147196E+00
1.122260E+00
1.106380E+00
1.124693E+00
@ -408,11 +403,6 @@ k cmfd
1.170308E+00
1.184540E+00
cmfd entropy
3.212002E+00
3.206393E+00
3.223984E+00
3.222764E+00
3.232123E+00
3.242083E+00
3.246067E+00
3.238869E+00
@ -425,11 +415,6 @@ cmfd entropy
3.221485E+00
3.219108E+00
cmfd balance
8.26180E-03
4.27338E-03
2.22686E-03
1.93026E-03
1.96979E-03
2.13756E-03
2.01479E-03
1.74519E-03
@ -442,11 +427,6 @@ cmfd balance
1.24780E-03
1.15560E-03
cmfd dominance ratio
5.404E-01
5.406E-01
5.449E-01
5.473E-01
5.534E-01
5.623E-01
5.738E-01
5.611E-01
@ -459,11 +439,6 @@ cmfd dominance ratio
5.412E-01
5.383E-01
cmfd openmc source comparison
1.575499E-02
1.293688E-02
3.531746E-03
8.281178E-03
5.771681E-03
7.459013E-03
5.012869E-03
1.770224E-03

View file

@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 10
cmfd_run.solver_begin = 10
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.window_type = 'expanding'

View file

@ -18,7 +18,7 @@ def test_cmfd_feed_ng():
cmfd_run.mesh = cmfd_mesh
cmfd_run.reset = [5]
cmfd_run.tally_begin = 10
cmfd_run.feedback_begin = 10
cmfd_run.solver_begin = 10
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.downscatter = True

View file

@ -391,11 +391,6 @@ cmfd indices
1.000000E+00
1.000000E+00
k cmfd
1.143785E+00
1.163460E+00
1.173453E+00
1.171056E+00
1.147214E+00
1.122230E+00
1.106385E+00
1.124706E+00
@ -408,11 +403,6 @@ k cmfd
1.170183E+00
1.184408E+00
cmfd entropy
3.211758E+00
3.206356E+00
3.223933E+00
3.222754E+00
3.232110E+00
3.242098E+00
3.246062E+00
3.238858E+00
@ -425,11 +415,6 @@ cmfd entropy
3.221498E+00
3.219126E+00
cmfd balance
8.26180E-03
4.27338E-03
2.22686E-03
1.93026E-03
1.96979E-03
2.13756E-03
2.01521E-03
1.74538E-03
@ -442,11 +427,6 @@ cmfd balance
1.24170E-03
1.15645E-03
cmfd dominance ratio
5.408E-01
5.374E-01
5.420E-01
5.444E-01
5.513E-01
5.613E-01
5.722E-01
5.592E-01
@ -459,11 +439,6 @@ cmfd dominance ratio
5.413E-01
5.381E-01
cmfd openmc source comparison
1.597982E-02
1.276493E-02
3.495229E-03
8.157807E-03
5.715330E-03
7.433111E-03
5.006211E-03
1.766072E-03

View file

@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 10
cmfd_run.solver_begin = 10
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.window_type = 'expanding'

View file

@ -391,11 +391,6 @@ cmfd indices
1.000000E+00
1.000000E+00
k cmfd
1.143597E+00
1.163387E+00
1.162391E+00
1.163351E+00
1.145721E+00
1.134785E+00
1.119048E+00
1.116124E+00
@ -408,11 +403,6 @@ k cmfd
1.166709E+00
1.167223E+00
cmfd entropy
3.212002E+00
3.206393E+00
3.222109E+00
3.221996E+00
3.229978E+00
3.234615E+00
3.246512E+00
3.244634E+00
@ -425,11 +415,6 @@ cmfd entropy
3.203758E+00
3.201798E+00
cmfd balance
8.26180E-03
4.27338E-03
2.62159E-03
2.40301E-03
2.08484E-03
1.58351E-03
1.59196E-03
1.87591E-03
@ -442,11 +427,6 @@ cmfd balance
2.25515E-03
1.53613E-03
cmfd dominance ratio
5.404E-01
5.406E-01
5.432E-01
5.454E-01
5.507E-01
5.578E-01
5.679E-01
5.671E-01
@ -459,11 +439,6 @@ cmfd dominance ratio
5.374E-01
5.321E-01
cmfd openmc source comparison
1.575499E-02
1.293688E-02
5.920734E-03
9.746850E-03
7.183896E-03
7.693485E-03
4.158805E-03
2.962505E-03

View file

@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 10
cmfd_run.solver_begin = 10
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.window_type = 'rolling'

View file

@ -15,7 +15,7 @@ def test_cmfd_nofeed():
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = False
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]

View file

@ -53,7 +53,7 @@ def test_cmfd_restart():
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.tally_begin = 5
cmfd_run.feedback_begin = 5
cmfd_run.solver_begin = 5
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run()
@ -62,7 +62,7 @@ def test_cmfd_restart():
cmfd_run2 = cmfd.CMFDRun()
cmfd_run2.mesh = cmfd_mesh2
cmfd_run2.tally_begin = 5
cmfd_run2.feedback_begin = 5
cmfd_run2.solver_begin = 5
cmfd_run2.feedback = True
cmfd_run2.gauss_seidel_tolerance = [1.e-15, 1.e-20]

View file

@ -63,4 +63,17 @@
<parameters>1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23</parameters>
</energy>
</source>
<source strength="0.1">
<space origin="1.0 1.0 0.0" type="cylindrical">
<r parameters="2.0 3.0" type="uniform" />
<phi parameters="0.0 6.283185307179586" type="uniform" />
<z interpolation="linear-linear" type="tabular">
<parameters>-2.0 0.0 2.0 0.2 0.3 0.2</parameters>
</z>
</space>
<angle type="isotropic" />
<energy interpolation="histogram" type="tabular">
<parameters>1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23</parameters>
</energy>
</source>
</settings>

View file

@ -1,2 +1,2 @@
k-combined:
3.004769E-01 3.944044E-03
2.999424E-01 9.520402E-03

View file

@ -38,7 +38,10 @@ class SourceTestHarness(PyAPITestHarness):
spatial3 = openmc.stats.Point([1.2, -2.3, 0.781])
spatial4 = openmc.stats.SphericalIndependent(r_dist, theta_dist,
phi_dist,
origin=(1.0, 1.0, 0.0))
origin=(1., 1., 0.))
spatial5 = openmc.stats.CylindricalIndependent(r_dist, phi_dist,
z_dist,
origin=(1., 1., 0.))
mu_dist = openmc.stats.Discrete([-1., 0., 1.], [0.5, 0.25, 0.25])
phi_dist = openmc.stats.Uniform(0., 6.28318530718)
@ -57,12 +60,13 @@ class SourceTestHarness(PyAPITestHarness):
source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3)
source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1)
source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1)
source5 = openmc.Source(spatial5, angle3, energy3, strength=0.1)
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 5
settings.particles = 1000
settings.source = [source1, source2, source3, source4]
settings.source = [source1, source2, source3, source4, source5]
settings.export_to_xml()

View file

@ -18,7 +18,7 @@ def test_data_library(tmpdir):
assert f['type'] == 'neutron'
assert 'U235' in f['materials']
f = lib.get_by_material('c_H_in_H2O')
f = lib.get_by_material('c_H_in_H2O', data_type='thermal')
assert f['type'] == 'thermal'
assert 'c_H_in_H2O' in f['materials']

View file

@ -7,7 +7,7 @@ if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
fi
# Download ENDF/B-VII.1 distribution
ENDF=$HOME/endf-b-vii.1/
ENDF=$HOME/endf-b-vii.1
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ
fi