'out_left': 0, 'in_left': 1, 'out_right': 2, 'in_right': 3,

'out_back': 4, 'in_back': 5, 'out_front': 6, 'in_front': 7,

'out_bottom': 8, 'in_bottom': 9, 'out_top': 10, 'in_top': 11

"""A structured Cartesian mesh used for CMFD acceleration.

def __init__(self):

self._lower_left = None

self._upper_right = None

self._dimension = None

self._width = None

self._energy = None

self._albedo = None

self._map = None

self._mesh_type = 'regular'

self._grid = None

def __repr__(self):

outstr = type(self).__name__ + '\n'

if self._mesh_type == 'regular':

outstr += (self._get_repr(self._lower_left, "Lower left") + "\n" +

self._get_repr(self._upper_right, "Upper right") + "\n" +

self._get_repr(self._dimension, "Dimension") + "\n" +

self._get_repr(self._width, "Width") + "\n" +

self._get_repr(self._albedo, "Albedo"))

elif self._mesh_type == 'rectilinear':

outstr += (self._get_repr(self._grid[0], "X-grid") + "\n" +

self._get_repr(self._grid[1], "Y-grid") + "\n" +

self._get_repr(self._grid[2], "Z-grid"))

return outstr

def _get_repr(self, list_var, label):

outstr = f"\t{label:<11} = "

if list(list_var):

outstr += ", ".join(str(i) for i in list_var)

return outstr

@property

def lower_left(self):

return self._lower_left

@property

def upper_right(self):

return self._upper_right

@property

def dimension(self):

return self._dimension

@property

def width(self):

return self._width

@property

def energy(self):

return self._energy

@property

def albedo(self):

return self._albedo

@property

def map(self):

return self._map

@property

def mesh_type(self):

return self._mesh_type

@property

def grid(self):

return self._grid

@lower_left.setter

def lower_left(self, lower_left):

check_type('CMFD mesh lower_left', lower_left, Iterable, Real)

check_length('CMFD mesh lower_left', lower_left, 2, 3)

self._lower_left = lower_left

self._display_mesh_warning('regular', 'CMFD mesh lower_left')

@upper_right.setter

def upper_right(self, upper_right):

check_type('CMFD mesh upper_right', upper_right, Iterable, Real)

check_length('CMFD mesh upper_right', upper_right, 2, 3)

self._upper_right = upper_right

self._display_mesh_warning('regular', 'CMFD mesh upper_right')

@dimension.setter

def dimension(self, dimension):

check_type('CMFD mesh dimension', dimension, Iterable, Integral)

check_length('CMFD mesh dimension', dimension, 2, 3)

for d in dimension:

check_greater_than('CMFD mesh dimension', d, 0)

self._dimension = dimension

@width.setter

def width(self, width):

check_type('CMFD mesh width', width, Iterable, Real)

check_length('CMFD mesh width', width, 2, 3)

for w in width:

check_greater_than('CMFD mesh width', w, 0)

self._width = width

self._display_mesh_warning('regular', 'CMFD mesh width')

@energy.setter

def energy(self, energy):

check_type('CMFD mesh energy', energy, Iterable, Real)

for e in energy:

check_greater_than('CMFD mesh energy', e, 0, True)

self._energy = energy

@albedo.setter

def albedo(self, albedo):

check_type('CMFD mesh albedo', albedo, Iterable, Real)

check_length('CMFD mesh albedo', albedo, 6)

for a in albedo:

check_greater_than('CMFD mesh albedo', a, 0, True)

check_less_than('CMFD mesh albedo', a, 1, True)

self._albedo = albedo

@map.setter

def map(self, mesh_map):

check_type('CMFD mesh map', mesh_map, Iterable, Integral)

for m in mesh_map:

check_value('CMFD mesh map', m, [0, 1])

self._map = mesh_map

@mesh_type.setter

def mesh_type(self, mesh_type):

check_value('CMFD mesh type', mesh_type, ['regular', 'rectilinear'])

self._mesh_type = mesh_type

@grid.setter

def grid(self, grid):

grid_length = 3

dims = ['x', 'y', 'z']

check_length('CMFD mesh grid', grid, grid_length)

for i in range(grid_length):

check_type(f'CMFD mesh {dims[i]}-grid', grid[i], Iterable,

Real)

check_greater_than(f'CMFD mesh {dims[i]}-grid length',

len(grid[i]), 1)

self._grid = [np.array(g) for g in grid]

self._display_mesh_warning('rectilinear', 'CMFD mesh grid')

def _display_mesh_warning(self, mesh_type, variable_label):

if self._mesh_type != mesh_type:

warn_msg = (f'Setting {variable_label} if mesh type is not set to '

f'{mesh_type} will have no effect')

r"""Class for running CMFD acceleration through the C API.

def __init__(self):

"""Constructor for CMFDRun class. Default values for instance variables

# Variables that users can modify

self._tally_begin = 1

self._solver_begin = 1

self._ref_d = np.array([])

self._display = {'balance': False, 'dominance': False,

'entropy': False, 'source': False}

self._downscatter = False

self._feedback = False

self._cmfd_ktol = 1.e-8

self._mesh = None

self._norm = 1.

self._power_monitor = False

self._run_adjoint = False

self._w_shift = 1.e6

self._stol = 1.e-8

self._reset = []

self._write_matrices = False

self._spectral = 0.0

self._gauss_seidel_tolerance = [1.e-10, 1.e-5]

self._adjoint_type = 'physical'

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

self._indices = np.zeros(4, dtype=np.int32)

self._egrid = None

self._albedo = None

self._coremap = None

self._mesh_id = None

self._tally_ids = None

self._energy_filters = None

self._cmfd_on = False

self._mat_dim = _CMFD_NOACCEL

self._keff_bal = None

self._keff = None

self._adj_keff = None

self._phi = None

self._adj_phi = None

self._openmc_src_rate = None

self._flux_rate = None

self._total_rate = None

self._p1scatt_rate = None

self._scatt_rate = None

self._nfiss_rate = None

self._current_rate = None

self._flux = None

self._totalxs = None

self._p1scattxs = None

self._scattxs = None

self._nfissxs = None

self._diffcof = None

self._dtilde = None

self._dhat = None

self._hxyz = None

self._current = None

self._cmfd_src = None

self._openmc_src = None

self._entropy = []

self._balance = []

self._src_cmp = []

self._dom = []

self._k_cmfd = []

self._resnb = None

self._reset_every = None

self._time_cmfd = None

self._time_cmfdbuild = None

self._time_cmfdsolve = None

# All index-related variables, for numpy vectorization

self._first_x_accel = None

self._last_x_accel = None

self._first_y_accel = None

self._last_y_accel = None

self._first_z_accel = None

self._last_z_accel = None

self._notfirst_x_accel = None

self._notlast_x_accel = None

self._notfirst_y_accel = None

self._notlast_y_accel = None

self._notfirst_z_accel = None

self._notlast_z_accel = None

self._is_adj_ref_left = None

self._is_adj_ref_right = None

self._is_adj_ref_back = None

self._is_adj_ref_front = None

self._is_adj_ref_bottom = None

self._is_adj_ref_top = None

self._accel_idxs = None

self._accel_neig_left_idxs = None

self._accel_neig_right_idxs = None

self._accel_neig_back_idxs = None

self._accel_neig_front_idxs = None

self._accel_neig_bot_idxs = None

self._accel_neig_top_idxs = None

self._loss_row = None

self._loss_col = None

self._prod_row = None

self._prod_col = None

@property

def tally_begin(self):

return self._tally_begin

@property

def solver_begin(self):

return self._solver_begin

@property

def ref_d(self):

return self._ref_d

@property

def display(self):

return self._display

@property

def downscatter(self):

return self._downscatter

@property

def feedback(self):

return self._feedback

@property

def cmfd_ktol(self):

return self._cmfd_ktol

@property

def mesh(self):

return self._mesh

@property

def norm(self):

return self._norm

@property

def adjoint_type(self):

return self._adjoint_type

@property

def window_type(self):

return self._window_type

@property

def window_size(self):

return self._window_size

@property

def power_monitor(self):

return self._power_monitor

@property

def run_adjoint(self):

return self._run_adjoint

@property

def w_shift(self):

return self._w_shift

@property

def stol(self):

return self._stol

@property

def spectral(self):

return self._spectral

@property

def reset(self):

return self._reset

@property

def write_matrices(self):

return self._write_matrices

@property

def gauss_seidel_tolerance(self):

return self._gauss_seidel_tolerance

@property

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

@property

def dom(self):

return self._dom

@property

def src_cmp(self):

return self._src_cmp

@property

def balance(self):

return self._balance

@property

def entropy(self):

return self._entropy

@property

def k_cmfd(self):

return self._k_cmfd

@tally_begin.setter

def tally_begin(self, begin):

check_type('CMFD tally begin batch', begin, Integral)

check_greater_than('CMFD tally begin batch', begin, 0)

self._tally_begin = 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._solver_begin = begin

@ref_d.setter

def ref_d(self, diff_params):

check_type('Reference diffusion params', diff_params,

Iterable, Real)

self._ref_d = np.array(diff_params)

@display.setter

def display(self, display):

check_type('display', display, Mapping)

for key, value in display.items():

check_value('display key', key,

('balance', 'entropy', 'dominance', 'source'))

check_type(f"display['{key}']", value, bool)

self._display[key] = value

@downscatter.setter

def downscatter(self, downscatter):

check_type('CMFD downscatter', downscatter, bool)

self._downscatter = downscatter

@feedback.setter

def feedback(self, feedback):

check_type('CMFD feedback', feedback, bool)

self._feedback = feedback

@cmfd_ktol.setter

def cmfd_ktol(self, cmfd_ktol):

check_type('CMFD eigenvalue tolerance', cmfd_ktol, Real)

self._cmfd_ktol = cmfd_ktol

@mesh.setter

def mesh(self, cmfd_mesh):

check_type('CMFD mesh', cmfd_mesh, CMFDMesh)

if cmfd_mesh.mesh_type == 'regular':

# Check dimension defined

if cmfd_mesh.dimension is None:

raise ValueError('CMFD regular mesh requires spatial '

'dimensions to be specified')

# Check lower left defined

if cmfd_mesh.lower_left is None:

raise ValueError('CMFD regular mesh requires lower left '

'coordinates to be specified')

# Check that both upper right and width both not defined

if cmfd_mesh.upper_right is not None and cmfd_mesh.width is not None:

raise ValueError('Both upper right coordinates and width '

'cannot be specified for CMFD regular mesh')

# Check that at least one of width or upper right is defined

if cmfd_mesh.upper_right is None and cmfd_mesh.width is None:

raise ValueError('CMFD regular mesh requires either upper right '

'coordinates or width to be specified')

# Check width and lower length are same dimension and define

# upper_right

if cmfd_mesh.width is not None:

check_length('CMFD mesh width', cmfd_mesh.width,

len(cmfd_mesh.lower_left))

cmfd_mesh.upper_right = np.array(cmfd_mesh.lower_left) + \

np.array(cmfd_mesh.width) * np.array(cmfd_mesh.dimension)

# Check upper_right and lower length are same dimension and define

# width

elif cmfd_mesh.upper_right is not None:

check_length('CMFD mesh upper right', cmfd_mesh.upper_right,

len(cmfd_mesh.lower_left))

# Check upper right coordinates are greater than lower left

if np.any(np.array(cmfd_mesh.upper_right) <=

np.array(cmfd_mesh.lower_left)):

raise ValueError('CMFD regular mesh requires upper right '

'coordinates to be greater than lower '

'left coordinates')

cmfd_mesh.width = np.true_divide((np.array(cmfd_mesh.upper_right) -

np.array(cmfd_mesh.lower_left)),

np.array(cmfd_mesh.dimension))

elif cmfd_mesh.mesh_type == 'rectilinear':

# Check dimension defined

if cmfd_mesh.grid is None:

raise ValueError('CMFD rectilinear mesh requires spatial '

'grid to be specified')

cmfd_mesh.dimension = [len(cmfd_mesh.grid[i]) - 1 for i in range(3)]

self._mesh = cmfd_mesh

@norm.setter

def norm(self, norm):

check_type('CMFD norm', norm, Real)

self._norm = norm

@adjoint_type.setter

def adjoint_type(self, adjoint_type):

check_type('CMFD adjoint type', adjoint_type, str)

check_value('CMFD adjoint type', adjoint_type,

['math', 'physical'])

self._adjoint_type = adjoint_type

@window_type.setter

def window_type(self, window_type):

check_type('CMFD window type', window_type, str)

check_value('CMFD window type', window_type,

['none', 'rolling', 'expanding'])

self._window_type = window_type

@window_size.setter

def window_size(self, window_size):

check_type('CMFD window size', window_size, Integral)

check_greater_than('CMFD window size', window_size, 0)

if self._window_type != 'rolling':

warn_msg = 'Window size will have no effect on CMFD simulation ' \

'unless window type is set to "rolling".'

warnings.warn(warn_msg, RuntimeWarning)

self._window_size = window_size

@power_monitor.setter

def power_monitor(self, power_monitor):

check_type('CMFD power monitor', power_monitor, bool)

self._power_monitor = power_monitor

@run_adjoint.setter

def run_adjoint(self, run_adjoint):

check_type('CMFD run adjoint', run_adjoint, bool)

self._run_adjoint = run_adjoint

@w_shift.setter

def w_shift(self, w_shift):

check_type('CMFD Wielandt shift', w_shift, Real)

self._w_shift = w_shift

@stol.setter

def stol(self, stol):

check_type('CMFD fission source tolerance', stol, Real)

self._stol = stol

@spectral.setter

def spectral(self, spectral):

check_type('CMFD spectral radius', spectral, Real)

self._spectral = spectral

@reset.setter

def reset(self, reset):

check_type('tally reset batches', reset, Iterable, Integral)

self._reset = reset

@write_matrices.setter

def write_matrices(self, write_matrices):

check_type('CMFD write matrices', write_matrices, bool)

self._write_matrices = write_matrices

@gauss_seidel_tolerance.setter

def gauss_seidel_tolerance(self, gauss_seidel_tolerance):

check_type('CMFD Gauss-Seidel tolerance', gauss_seidel_tolerance,

Iterable, Real)

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

with self.run_in_memory(**kwargs):

for _ in self.iter_batches():

pass

@contextmanager

def run_in_memory(self, **kwargs):

""" Context manager for running CMFD functions with OpenMC shared

# Store intracomm for part of CMFD routine where MPI reduce and

# broadcast calls are made

if 'intracomm' in kwargs and kwargs['intracomm'] is not None:

self._intracomm = kwargs['intracomm']

elif have_mpi:

self._intracomm = MPI.COMM_WORLD

# Run and pass arguments to C API run_in_memory function

with openmc.lib.run_in_memory(**kwargs):

self.init()

yield

self.finalize()

def iter_batches(self):

""" Iterator over batches.

status = 0

while status == 0:

status = self.next_batch()

yield

def init(self):

""" Initialize CMFDRun instance by setting up CMFD parameters and

# Configure CMFD parameters

self._configure_cmfd()

# Create tally objects

self._create_cmfd_tally()

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

# Initialize all variables used for linear solver in C++

self._initialize_linsolver()

# Initialize simulation

openmc.lib.simulation_init()

# Set cmfd_run variable to True through C API

openmc.lib.settings.cmfd_run = True

def next_batch(self):

""" Run next batch for CMFDRun.

# Initialize CMFD batch

self._cmfd_init_batch()

# Run next batch

status = openmc.lib.next_batch()

# Perform CMFD calculations

self._execute_cmfd()

# Write CMFD data to statepoint

if openmc.lib.is_statepoint_batch():

self.statepoint_write()

return status

def finalize(self):

""" Finalize simulation by calling

# Finalize simulation

openmc.lib.simulation_finalize()

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

if filename is None:

batch_str_len = len(str(openmc.lib.settings.get_batches()))

batch_str = str(openmc.lib.current_batch()).zfill(batch_str_len)

filename = f'statepoint.{batch_str}.h5'

# Call C API statepoint_write to save source distribution with CMFD

# feedback

openmc.lib.statepoint_write(filename=filename)

# Append CMFD data to statepoint file using h5py

self._write_cmfd_statepoint(filename)

def _write_cmfd_statepoint(self, filename):

"""Append all CNFD simulation parameters to existing statepoint

if openmc.lib.master():

with h5py.File(filename, 'a') as f:

if 'cmfd' not in f:

if openmc.lib.settings.verbosity >= 5:

print(f' Writing CMFD data to {filename}...')

sys.stdout.flush()

cmfd_group = f.create_group("cmfd")

cmfd_group.attrs['cmfd_on'] = self._cmfd_on

cmfd_group.attrs['feedback'] = self._feedback

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

cmfd_group.attrs['time_cmfdbuild'] = self._time_cmfdbuild

cmfd_group.attrs['time_cmfdsolve'] = self._time_cmfdsolve

cmfd_group.attrs['window_size'] = self._window_size

cmfd_group.attrs['window_type'] = self._window_type

cmfd_group.create_dataset('k_cmfd', data=self._k_cmfd)

cmfd_group.create_dataset('dom', data=self._dom)

cmfd_group.create_dataset('src_cmp', data=self._src_cmp)

cmfd_group.create_dataset('balance', data=self._balance)

cmfd_group.create_dataset('entropy', data=self._entropy)

cmfd_group.create_dataset('reset', data=self._reset)

cmfd_group.create_dataset('albedo', data=self._albedo)

cmfd_group.create_dataset('coremap', data=self._coremap)

cmfd_group.create_dataset('egrid', data=self._egrid)

cmfd_group.create_dataset('indices', data=self._indices)

cmfd_group.create_dataset('tally_ids',

data=self._tally_ids)

cmfd_group.create_dataset('current_rate',

data=self._current_rate)

cmfd_group.create_dataset('flux_rate',

data=self._flux_rate)

cmfd_group.create_dataset('nfiss_rate',

data=self._nfiss_rate)

cmfd_group.create_dataset('openmc_src_rate',

data=self._openmc_src_rate)

cmfd_group.create_dataset('p1scatt_rate',

data=self._p1scatt_rate)

cmfd_group.create_dataset('scatt_rate',

data=self._scatt_rate)

cmfd_group.create_dataset('total_rate',

data=self._total_rate)

elif openmc.settings.verbosity >= 5:

print(' CMFD data not written to statepoint file as it '

'already exists in {}'.format(filename), flush=True)

def _initialize_linsolver(self):

# Determine number of rows in CMFD matrix

ng = self._indices[3]

n = self._mat_dim*ng

# Create temp loss matrix to pass row/col indices to C++ linear solver

loss_row = self._loss_row

loss_col = self._loss_col

temp_data = np.ones(len(loss_row))

temp_loss = csr_array((temp_data, (loss_row, loss_col)), shape=(n, n))

temp_loss.sort_indices()

# Pass coremap as 1-d array of 32-bit integers

coremap = np.swapaxes(self._coremap, 0, 2).flatten().astype(np.int32)

return openmc.lib._dll.openmc_initialize_linsolver(

temp_loss.indptr.astype(np.int32), len(temp_loss.indptr),

temp_loss.indices.astype(np.int32), len(temp_loss.indices), n,

self._spectral, coremap, self._use_all_threads

)

def _write_cmfd_output(self):

"""Write CMFD output to buffer at the end of each batch"""

# Display CMFD k-effective

outstr = '{:>11s}CMFD k: {:0.5f}'.format('', self._k_cmfd[-1])

# Display value of additional fields based on display dict