"""This module can be used to specify parameters used for coarse mesh finite
difference (CMFD) acceleration in OpenMC. CMFD was first proposed by [Smith]_
and is widely used in accelerating neutron transport problems.

References
----------

.. [Smith] K. Smith, "Nodal method storage reduction by non-linear
   iteration", *Trans. Am. Nucl. Soc.*, **44**, 265 (1983).

"""

from collections.abc import Iterable, Mapping
from contextlib import contextmanager
from numbers import Real, Integral
import sys
import time
import warnings

import h5py
import numpy as np
from scipy import sparse

import openmc.lib
from .checkvalue import (check_type, check_length, check_value,
                         check_greater_than, check_less_than)
from .exceptions import OpenMCError

# See if mpi4py module can be imported, define have_mpi global variable
try:
    from mpi4py import MPI
    have_mpi = True
except ImportError:
    have_mpi = False

# Maximum/minimum neutron energies
_ENERGY_MAX_NEUTRON = np.inf
_ENERGY_MIN_NEUTRON = 0.

# Tolerance for detecting zero flux values
_TINY_BIT = 1.e-8

# For non-accelerated regions on coarse mesh overlay
_CMFD_NOACCEL = -1

# Constant to represent a zero flux "albedo"
_ZERO_FLUX = 999.0

# Map that returns index of current direction in numpy current matrix
_CURRENTS = {
    '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
}


class CMFDMesh:
    """A structured Cartesian mesh used for CMFD acceleration.

    Attributes
    ----------
    lower_left : Iterable of float
        The lower-left corner of a regular structured mesh. If only two
        coordinates are given, it is assumed that the mesh is an x-y mesh.
    upper_right : Iterable of float
        The upper-right corner of a regular structured mesh. If only two
        coordinates are given, it is assumed that the mesh is an x-y mesh.
    dimension : Iterable of int
        The number of mesh cells in each direction for a regular structured
        mesh.
    width : Iterable of float
        The width of mesh cells in each direction for a regular structured
        mesh
    energy : Iterable of float
        Energy bins in eV, listed in ascending order (e.g. [0.0, 0.625e-1,
        20.0e6]) for CMFD tallies and acceleration. If no energy bins are
        listed, OpenMC automatically assumes a one energy group calculation
        over the entire energy range.
    albedo : Iterable of float
        Surface ratio of incoming to outgoing partial currents on global
        boundary conditions. They are listed in the following order: -x +x -y
        +y -z +z.
    map : Iterable of int
        An optional acceleration map can be specified to overlay on the coarse
        mesh spatial grid. If this option is used, a ``0`` is used for a
        non-accelerated region and a ``1`` is used for an accelerated region.
        For a simple 4x4 coarse mesh with a 2x2 fuel lattice surrounded by
        reflector, the map is:

        ::

            [0, 0, 0, 0,
             0, 1, 1, 0,
             0, 1, 1, 0,
             0, 0, 0, 0]

        Therefore a 2x2 system of equations is solved rather than a 4x4. This
        is extremely important to use in reflectors as neutrons will not
        contribute to any tallies far away from fission source neutron regions.
        A ``1`` must be used to identify any fission source region.
    mesh_type : str
        Type of structured mesh to use. Acceptable values are:
        * "regular" - Use RegularMesh to define CMFD mesh
        * "rectilinear" - Use RectilinearMesh to define CMFD
    grid : Iterable of Iterable of float
        Grid used to define RectilinearMesh. First dimension must have length
        3 where grid[0], grid[1], and grid[2] correspond to the x-, y-, and
        z-grids respectively

    """

    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 = "\t{:<11} = ".format(label)
        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('CMFD mesh {}-grid'.format(dims[i]), grid[i], Iterable,
                       Real)
            check_greater_than('CMFD mesh {}-grid length'.format(dims[i]),
                               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')
            warnings.warn(warn_msg, RuntimeWarning)


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

    Attributes
    ----------
    tally_begin : int
        Batch number at which CMFD tallies should begin accumulating
    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
    display : dict
        Dictionary indicating which CMFD results to output. Note that CMFD
        k-effective will always be outputted. Acceptable keys are:

        * "balance" - Whether to output RMS [%] of the residual from the
          neutron balance equation on CMFD tallies (bool)
        * "dominance" - Whether to output the estimated dominance ratio from
          the CMFD iterations (bool)
        * "entropy" - Whether to output the *entropy* of the CMFD predicted
          fission source (bool)
        * "source" - Whether to output the RMS [%] between the OpenMC fission
          source and CMFD fission source (bool)

    downscatter : bool
        Indicate whether an effective downscatter cross section should be used
        when using 2-group CMFD.
    feedback : bool
        Indicate whether or not the CMFD diffusion result is used to adjust the weight
        of fission source neutrons on the next OpenMC batch. Defaults to False.
    cmfd_ktol : float
        Tolerance on the eigenvalue when performing CMFD power iteration
    mesh : openmc.cmfd.CMFDMesh
        Structured mesh to be used for acceleration
    norm : float
        Normalization factor applied to the CMFD fission source distribution
    power_monitor : bool
        View convergence of power iteration during CMFD acceleration
    run_adjoint : bool
        Perform adjoint calculation on the last batch
    w_shift : float
        Optional Wielandt shift parameter for accelerating power iterations. By
        default, it is very large so there is effectively no impact.
    stol : float
        Tolerance on the fission source when performing CMFD power iteration
    reset : list of int
        List of batch numbers at which CMFD tallies should be reset
    write_matrices : bool
        Write sparse matrices that are used during CMFD acceleration (loss,
        production) and resultant normalized flux vector phi to file
    spectral : float
        Optional spectral radius that can be used to accelerate the convergence
        of Gauss-Seidel iterations during CMFD power iteration.
    gauss_seidel_tolerance : Iterable of float
        Two parameters specifying the absolute inner tolerance and the relative
        inner tolerance for Gauss-Seidel iterations when performing CMFD.
    adjoint_type : {'physical', 'math'}
        Stores type of adjoint calculation that should be performed.
        ``run_adjoint`` must be true for an adjoint calculation to be
        performed. Options are:

        * "physical" - Create adjoint matrices from physical parameters of
          CMFD problem
        * "math" - Create adjoint matrices mathematically as the transpose of
          loss and production CMFD matrices

    window_type : {'expanding', 'rolling', 'none'}
        Specifies type of tally window scheme to use to accumulate CMFD
        tallies. Options are:

          * "expanding" - Have an expanding window that doubles in size
            to give more weight to more recent tallies as more generations are
            simulated
          * "rolling" - Have a fixed window size that aggregates tallies from
            the same number of previous generations tallied
          * "none" - Don't use a windowing scheme so that all tallies from last
            time they were reset are used for the CMFD algorithm.

    window_size : int
        Size of window to use for tally window scheme. Only relevant when
        window_type is set to "rolling"
    indices : numpy.ndarray
        Stores spatial and group dimensions as [nx, ny, nz, ng]
    cmfd_src : numpy.ndarray
        CMFD source distribution calculated from solving CMFD equations
    entropy : list of floats
        "Shannon entropy" from CMFD fission source, stored for each generation
        that CMFD is invoked
    balance : list of floats
        RMS of neutron balance equations, stored for each generation that CMFD
        is invoked
    src_cmp : list of floats
        RMS deviation of OpenMC and CMFD normalized source, stored for each
        generation that CMFD is invoked
    dom : list of floats
        Dominance ratio from solving CMFD matrix equations, stored for each
        generation that CMFD is invoked
    k_cmfd : list of floats
        List of CMFD k-effectives, stored for each generation that CMFD is
        invoked
    time_cmfd : float
        Time for entire CMFD calculation, in seconds
    time_cmfdbuild : float
        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

    """

    def __init__(self):
        """Constructor for CMFDRun class. Default values for instance variables
        set in this method.

        """
        # 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("display['{}']".format(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

        This method is called by the user to run CMFD once instance variables of
        CMFDRun class are set

        Parameters
        ----------
        **kwargs
            All keyword arguments are passed to
            :func:`openmc.lib.run_in_memory`.

        """
        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
        library functions.

        This function can be used with a 'with' statement to ensure the
        CMFDRun class is properly initialized/finalized. For example::

            from openmc import cmfd
            cmfd_run = cmfd.CMFDRun()
            with cmfd_run.run_in_memory():
                do_stuff_before_simulation_start()
                for _ in cmfd_run.iter_batches():
                    do_stuff_between_batches()

        Parameters
        ----------
        **kwargs
            All keyword arguments passed to :func:`openmc.lib.run_in_memory`.

        """
        # 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.

        This function returns a generator-iterator that allows Python code to
        be run between batches when running an OpenMC simulation with CMFD.
        It should be used in conjunction with
        :func`openmc.cmfd.CMFDRun.run_in_memory` to ensure proper
        initialization/finalization of CMFDRun instance.

        """
        status = 0
        while status == 0:
            status = self.next_batch()
            yield

    def init(self):
        """ Initialize CMFDRun instance by setting up CMFD parameters and
        calling :func:`openmc.lib.simulation_init`

        """
        # 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.

        Returns
        -------
        int
            Status after running a batch (0=normal, 1=reached maximum number of
            batches, 2=tally triggers reached)

        """
        # 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
        :func:`openmc.lib.simulation_finalize` and print out CMFD timing
        information.

        """
        # 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

        Parameters
        ----------
        filename : str
            Filename of 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

        Parameters
        ----------
        filename : str
            Filename of statepoint

        """
        if openmc.lib.master():
            with h5py.File(filename, 'a') as f:
                if 'cmfd' not in f:
                    if openmc.lib.settings.verbosity >= 5:
                        print(' Writing CMFD data to {}...'.format(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 = sparse.csr_matrix((temp_data, (loss_row, loss_col)),
                                      shape=(n, n))

        # Pass coremap as 1-d array of 32-bit integers
        coremap = np.swapaxes(self._coremap, 0, 2).flatten().astype(np.int32)

        args = temp_loss.indptr, len(temp_loss.indptr), \
            temp_loss.indices, len(temp_loss.indices), n, \
            self._spectral, coremap, self._use_all_threads
        return openmc.lib._dll.openmc_initialize_linsolver(*args)

    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
        outstr += '\n'
        if self._display['dominance']:
            outstr += ('{:>11s}Dom Rat:   {:0.5f}\n'
                       .format('', self._dom[-1]))
        if self._display['entropy']:
            outstr += ('{:>11s}CMFD Ent:  {:0.5f}\n'
                       .format('', self._entropy[-1]))
        if self._display['source']:
            outstr += ('{:>11s}RMS Src:   {:0.5f}\n'
                       .format('', self._src_cmp[-1]))
        if self._display['balance']:
            outstr += ('{:>11s}RMS Bal:   {:0.5f}\n'
                       .format('', self._balance[-1]))

        print(outstr)
        sys.stdout.flush()

    def _write_cmfd_timing_stats(self):
        """Write CMFD timing stats to buffer after finalizing simulation"""
        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:
            # 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 _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')
            sys.stdout.flush()

        # Check if CMFD mesh is defined
        if self._mesh is None:
            raise ValueError('No CMFD mesh has been specified for '
                             'simulation')

        # Set spatial dimensions of CMFD object
        for i, n in enumerate(self._mesh.dimension):
            self._indices[i] = n

        # Check if in continuous energy mode
        if not openmc.lib.settings.run_CE:
            raise OpenMCError('CMFD must be run in continuous energy mode')

        # Set number of energy groups
        if self._mesh.energy is not None:
            ng = len(self._mesh.energy)
            self._egrid = np.array(self._mesh.energy)
            self._indices[3] = ng - 1
            self._energy_filters = True
        else:
            self._egrid = np.array([_ENERGY_MIN_NEUTRON, _ENERGY_MAX_NEUTRON])
            self._indices[3] = 1
            self._energy_filters = False

        # 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.prod(self._indices[:3]))
            if openmc.lib.master():
                self._coremap = np.array(self._mesh.map)
        else:
            if openmc.lib.master():
                self._coremap = np.ones(np.prod(self._indices[:3]), dtype=int)

        # Check CMFD tallies accummulated before feedback turned on
        if self._feedback and self._solver_begin < self._tally_begin:
            raise ValueError('Tally begin must be less than or equal to '
                             'CMFD begin')

        # Initialize parameters for CMFD tally windows
        self._set_tally_window()

        # Extract spatial and energy indices
        nx, ny, nz, ng = self._indices

        # Initialize CMFD source to all ones
        self._cmfd_src = np.ones((nx, ny, nz, ng))

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

            # 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

        Parameters
        ----------
        filename : str
            Filename of statepoint to read from

        """
        with h5py.File(filename, 'r') as f:
            if 'cmfd' not in f:
                raise OpenMCError('Could not find CMFD parameters in ',
                                  'file {}'.format(filename))
            else:
                # Overwrite CMFD values from statepoint
                if (openmc.lib.master() and
                        openmc.lib.settings.verbosity >= 5):
                    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._solver_begin = cmfd_group.attrs['solver_begin']
                self._tally_begin = cmfd_group.attrs['tally_begin']
                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._egrid = cmfd_group['egrid'][()]
                self._indices = cmfd_group['indices'][()]
                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'])
                cmfd_mesh = f['tallies']['meshes'][cmfd_mesh_name]
                self._mesh.mesh_type = cmfd_mesh['type'][()].decode()
                if self._mesh.mesh_type == 'regular':
                    self._mesh.dimension = cmfd_mesh['dimension'][()]
                    self._mesh.lower_left = cmfd_mesh['lower_left'][()]
                    self._mesh.upper_right = cmfd_mesh['upper_right'][()]
                    self._mesh.width = cmfd_mesh['width'][()]
                elif self._mesh.mesh_type == 'rectilinear':
                    x_grid = cmfd_mesh['x_grid'][()]
                    y_grid = cmfd_mesh['y_grid'][()]
                    z_grid = cmfd_mesh['z_grid'][()]
                    self._mesh.grid = [x_grid, y_grid, z_grid]

                # 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"""
        # Set parameters for expanding window
        if self._window_type == 'expanding':
            self._reset_every = True
            self._window_size = 1
        # Set parameters for rolling window
        elif self.window_type == 'rolling':
            self._reset_every = True
        # Set parameters for default case, with no window
        else:
            self._window_size = 1
            self._reset_every = False

    def _cmfd_init_batch(self):
        """Handles CMFD options at the beginning of each batch"""
        # Get current batch through C API
        # Add 1 as next_batch has not been called yet
        current_batch = openmc.lib.current_batch() + 1

        # Check to activate CMFD solver and possible feedback
        if self._solver_begin == current_batch:
            self._cmfd_on = True

        # Check to reset tallies
        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"""
        if openmc.lib.master():
            # Start CMFD timer
            time_start_cmfd = time.time()

            if openmc.lib.current_batch() >= self._tally_begin:
                # Calculate all cross sections based on tally window averages
                self._compute_xs()

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

                # Call solver
                self._cmfd_solver_execute()

                # Store k-effective
                self._k_cmfd.append(self._keff)

                # Check to perform adjoint on last batch
                batches = openmc.lib.settings.get_batches()
                if openmc.lib.current_batch() == batches and self._run_adjoint:
                    self._cmfd_solver_execute(adjoint=True)

                # Calculate fission source
                self._calc_fission_source()

            # 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():
            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"""
        # Print message
        if (openmc.lib.settings.verbosity >= 6 and openmc.lib.master() and
                not self._reset_every):
            print(' CMFD tallies reset')
            sys.stdout.flush()

        # Reset CMFD tallies
        tallies = openmc.lib.tallies
        for tally_id in self._tally_ids:
            tallies[tally_id].reset()

    def _set_up_cmfd(self):
        """Configures CMFD object for a CMFD eigenvalue calculation

        """
        # Compute effective downscatter cross section
        if self._downscatter:
            self._compute_effective_downscatter()

        # Check neutron balance
        self._neutron_balance()

        # Calculate dtilde
        self._compute_dtilde()

        # Calculate dhat
        self._compute_dhat()

    def _cmfd_solver_execute(self, adjoint=False):
        """Sets up and runs power iteration solver for CMFD

        Parameters
        ----------
        adjoint : bool
            Whether or not to run an adjoint calculation

        """
        # Start timer for build
        time_start_buildcmfd = time.time()

        # Build the loss and production matrices
        if not adjoint:
            # Build matrices without adjoint calculation
            loss = self._build_loss_matrix(False)
            prod = self._build_prod_matrix(False)
        else:
            # Build adjoint matrices by running adjoint calculation
            if self._adjoint_type == 'physical':
                loss = self._build_loss_matrix(True)
                prod = self._build_prod_matrix(True)
            # Build adjoint matrices as transpose of non-adjoint matrices
            else:
                loss = self._build_loss_matrix(False).transpose()
                prod = self._build_prod_matrix(False).transpose()

        # Write out the matrices.
        if self._write_matrices:
            if not adjoint:
                self._write_matrix(loss, 'loss')
                self._write_matrix(prod, 'prod')
            else:
                self._write_matrix(loss, 'adj_loss')
                self._write_matrix(prod, 'adj_prod')

        # Stop timer for build
        time_stop_buildcmfd = time.time()
        self._time_cmfdbuild += time_stop_buildcmfd - time_start_buildcmfd

        # Begin power iteration
        time_start_solvecmfd = time.time()
        phi, keff, dom = self._execute_power_iter(loss, prod)
        time_stop_solvecmfd = time.time()
        self._time_cmfdsolve += time_stop_solvecmfd - time_start_solvecmfd

        # Save results, normalizing phi to sum to 1
        if adjoint:
            self._adj_keff = keff
            self._adj_phi = phi/np.sqrt(np.sum(phi*phi))
        else:
            self._keff = keff
            self._phi = phi/np.sqrt(np.sum(phi*phi))

        self._dom.append(dom)

        # Write out flux vector
        if self._write_matrices:
            if adjoint:
                self._write_vector(self._adj_phi, 'adj_fluxvec')
            else:
                self._write_vector(self._phi, 'fluxvec')

    def _write_vector(self, vector, base_filename):
        """Write a 1-D numpy array to file and also save it in .npy format.
        This particular format allows users to load the variable directly in a
        Python session with np.load()

        Parameters
        ----------
        vector : numpy.ndarray
            Vector that will be saved
        base_filename : str
            Filename to save vector as, without any file extension at the end.
            Vector will be saved to file [base_filename].dat and in numpy
            format as [base_filename].npy

        """
        # Write each element in vector to file
        np.savetxt(f'{base_filename}.dat', vector, fmt='%.8f')

        # Save as numpy format
        np.save(base_filename, vector)

    def _write_matrix(self, matrix, base_filename):
        """Write a numpy matrix to file and also save it in .npz format. This
        particular format allows users to load the variable directly in a
        Python session with scipy.sparse.load_npz()

        Parameters
        ----------
        matrix : scipy.sparse.spmatrix
            Sparse matrix that will be saved
        base_filename : str
            Filename to save matrix entries, without any file extension at the
            end. Matrix entries will be saved to file [base_filename].dat and
            in scipy format as [base_filename].npz

        """
        # Write row, col, and data of each entry in sparse matrix. This ignores
        # all zero-entries, and indices are written with zero-based indexing
        with open(base_filename+'.dat', 'w') as fh:
            for row in range(matrix.shape[0]):
                # Get all cols for particular row in matrix
                cols = matrix.indices[matrix.indptr[row]:matrix.indptr[row+1]]
                # Get all data entries for particular row in matrix
                data = matrix.data[matrix.indptr[row]:matrix.indptr[row+1]]
                for i in range(len(cols)):
                    fh.write('{:3d}, {:3d}, {:0.8f}\n'.format(
                        row, cols[i], data[i]))

        # Save matrix in scipy format
        sparse.save_npz(base_filename, matrix)

    def _calc_fission_source(self):
        """Calculates CMFD fission source from CMFD flux. If a coremap is
        defined, there will be a discrepancy between the spatial indices in the
        variables ``phi`` and ``nfissxs``, so ``phi`` needs to be mapped to the
        spatial indices of the cross sections. This can be done in a vectorized
        numpy manner or with for loops

        """
        # Extract number of groups and number of accelerated regions
        nx, ny, nz, ng = self._indices
        n = self._mat_dim

        # Compute cmfd_src in a vectorized manner by phi to the spatial
        # indices of the actual problem so that cmfd_flux can be multiplied by
        # nfissxs

        # Reshape phi by number of groups
        phi = self._phi.reshape((n, ng))

        # Extract indices of coremap that are accelerated
        idx = self._accel_idxs

        # Initialize CMFD flux map that maps phi to actual spatial and
        # group indices of problem
        cmfd_flux = np.zeros((nx, ny, nz, ng))

        # Loop over all groups and set CMFD flux based on indices of
        # coremap and values of phi
        for g in range(ng):
            phi_g = phi[:,g]
            cmfd_flux[idx + (g,)] = phi_g[self._coremap[idx]]

        # Compute fission source
        cmfd_src = (np.sum(self._nfissxs[:,:,:,:,:] *
                    cmfd_flux[:,:,:,:,np.newaxis], axis=3))

        # Normalize source such that it sums to 1.0
        self._cmfd_src = cmfd_src / np.sum(cmfd_src)

        # Compute entropy
        if openmc.lib.settings.entropy_on:
            # Compute source times log_2(source)
            source = self._cmfd_src[self._cmfd_src > 0] \
                * np.log(self._cmfd_src[self._cmfd_src > 0])/np.log(2)

            # Sum source and store
            self._entropy.append(-1.0 * np.sum(source))

        # Normalize source so average is 1.0
        self._cmfd_src = self._cmfd_src/np.sum(self._cmfd_src) * self._norm

        # Calculate differences between normalized sources
        self._src_cmp.append(np.sqrt(1.0 / self._norm
                             * np.sum((self._cmfd_src - self._openmc_src)**2)))

    def _build_loss_matrix(self, adjoint):
        # Extract spatial and energy indices and define matrix dimension
        ng = self._indices[3]
        n = self._mat_dim*ng

        # Define data entries used to build csr matrix
        data = np.array([])

        dtilde_left = self._dtilde[:,:,:,:,0]
        dtilde_right = self._dtilde[:,:,:,:,1]
        dtilde_back = self._dtilde[:,:,:,:,2]
        dtilde_front = self._dtilde[:,:,:,:,3]
        dtilde_bottom = self._dtilde[:,:,:,:,4]
        dtilde_top = self._dtilde[:,:,:,:,5]
        dhat_left = self._dhat[:,:,:,:,0]
        dhat_right = self._dhat[:,:,:,:,1]
        dhat_back = self._dhat[:,:,:,:,2]
        dhat_front = self._dhat[:,:,:,:,3]
        dhat_bottom = self._dhat[:,:,:,:,4]
        dhat_top = self._dhat[:,:,:,:,5]

        dx = self._hxyz[:,:,:,np.newaxis,0]
        dy = self._hxyz[:,:,:,np.newaxis,1]
        dz = self._hxyz[:,:,:,np.newaxis,2]

        # Define net leakage coefficient for each surface in each matrix
        # element
        jnet = (((dtilde_right + dhat_right)-(-1.0 * dtilde_left + dhat_left))
                / dx +
                ((dtilde_front + dhat_front)-(-1.0 * dtilde_back + dhat_back))
                / dy +
                ((dtilde_top + dhat_top)-(-1.0 * dtilde_bottom + dhat_bottom))
                / dz)

        for g in range(ng):
            # Define leakage terms that relate terms to their neighbors to the
            # left
            dtilde = self._dtilde[:,:,:,g,0][self._accel_neig_left_idxs]
            dhat = self._dhat[:,:,:,g,0][self._accel_neig_left_idxs]
            dx = self._hxyz[:,:,:,0][self._accel_neig_left_idxs]
            vals = (-1.0 * dtilde - dhat) / dx
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define leakage terms that relate terms to their neighbors to the
            # right
            dtilde = self._dtilde[:,:,:,g,1][self._accel_neig_right_idxs]
            dhat = self._dhat[:,:,:,g,1][self._accel_neig_right_idxs]
            dx = self._hxyz[:,:,:,0][self._accel_neig_right_idxs]
            vals = (-1.0 * dtilde + dhat) / dx
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define leakage terms that relate terms to their neighbors in the
            # back
            dtilde = self._dtilde[:,:,:,g,2][self._accel_neig_back_idxs]
            dhat = self._dhat[:,:,:,g,2][self._accel_neig_back_idxs]
            dy = self._hxyz[:,:,:,1][self._accel_neig_back_idxs]
            vals = (-1.0 * dtilde - dhat) / dy
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define leakage terms that relate terms to their neighbors in the
            # front
            dtilde = self._dtilde[:,:,:,g,3][self._accel_neig_front_idxs]
            dhat = self._dhat[:,:,:,g,3][self._accel_neig_front_idxs]
            dy = self._hxyz[:,:,:,1][self._accel_neig_front_idxs]
            vals = (-1.0 * dtilde + dhat) / dy
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define leakage terms that relate terms to their neighbors to the
            # bottom
            dtilde = self._dtilde[:,:,:,g,4][self._accel_neig_bot_idxs]
            dhat = self._dhat[:,:,:,g,4][self._accel_neig_bot_idxs]
            dz = self._hxyz[:,:,:,2][self._accel_neig_bot_idxs]
            vals = (-1.0 * dtilde - dhat) / dz
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define leakage terms that relate terms to their neighbors to the
            # top
            dtilde = self._dtilde[:,:,:,g,5][self._accel_neig_top_idxs]
            dhat = self._dhat[:,:,:,g,5][self._accel_neig_top_idxs]
            dz = self._hxyz[:,:,:,2][self._accel_neig_top_idxs]
            vals = (-1.0 * dtilde + dhat) / dz
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define terms that relate to loss of neutrons in a cell. These
            # correspond to all the diagonal entries of the loss matrix
            jnet_g = jnet[:,:,:,g][self._accel_idxs]
            total_xs = self._totalxs[:,:,:,g][self._accel_idxs]
            scatt_xs = self._scattxs[:,:,:,g,g][self._accel_idxs]
            vals = jnet_g + total_xs - scatt_xs
            # Store data to add to CSR matrix
            data = np.append(data, vals)

            # Define terms that relate to in-scattering from group to group.
            # These terms relate a mesh index to all mesh indices with the same
            # spatial dimensions but belong to a different energy group
            for h in range(ng):
                if h != g:
                    # Get scattering macro xs, transposed
                    if adjoint:
                        scatt_xs = self._scattxs[:,:,:,g,h][self._accel_idxs]
                    # Get scattering macro xs
                    else:
                        scatt_xs = self._scattxs[:,:,:,h,g][self._accel_idxs]
                    vals = -1.0 * scatt_xs
                    # Store data to add to CSR matrix
                    data = np.append(data, vals)

        # Create csr matrix
        loss_row = self._loss_row
        loss_col = self._loss_col
        loss = sparse.csr_matrix((data, (loss_row, loss_col)), shape=(n, n))
        return loss

    def _build_prod_matrix(self, adjoint):
        # Extract spatial and energy indices and define matrix dimension
        ng = self._indices[3]
        n = self._mat_dim*ng

        # Define rows, columns, and data used to build csr matrix
        data = np.array([])

        # Define terms that relate to fission production from group to group.
        for g in range(ng):
            for h in range(ng):
                # Get nu-fission macro xs, transposed
                if adjoint:
                    vals = (self._nfissxs[:, :, :, g, h])[self._accel_idxs]
                # Get nu-fission macro xs
                else:
                    vals = (self._nfissxs[:, :, :, h, g])[self._accel_idxs]
                # Store rows, cols, and data to add to CSR matrix
                data = np.append(data, vals)

        # Create csr matrix
        prod_row = self._prod_row
        prod_col = self._prod_col
        prod = sparse.csr_matrix((data, (prod_row, prod_col)), shape=(n, n))
        return prod

    def _execute_power_iter(self, loss, prod):
        """Main power iteration routine for the CMFD calculation

        Parameters
        ----------
        loss : scipy.sparse.spmatrix
            Sparse matrix storing elements of CMFD loss matrix
        prod : scipy.sparse.spmatrix
            Sparse matrix storing elements of CMFD production matrix

        Returns
        -------
        phi_n : numpy.ndarray
            Flux vector of CMFD problem
        k_n : float
            Eigenvalue of CMFD problem
        dom : float
            Dominance ratio of CMFD problem

        """
        # Get problem size
        n = loss.shape[0]

        # Set up tolerances for C++ solver
        atoli = self._gauss_seidel_tolerance[0]
        rtoli = self._gauss_seidel_tolerance[1]
        toli = rtoli * 100

        # Set up flux vectors, intital guess set to 1
        phi_n = np.ones((n,))
        phi_o = np.ones((n,))

        # Set up source vectors
        s_n = np.zeros((n,))
        s_o = np.zeros((n,))

        # Set initial guess
        k_n = openmc.lib.keff()[0]
        k_o = k_n
        dw = self._w_shift
        k_s = k_o + dw
        k_ln = 1.0/(1.0/k_n - 1.0/k_s)
        k_lo = k_ln

        # Set norms to 0
        norm_n = 0.0
        norm_o = 0.0

        # Maximum number of power iterations
        maxits = 10000

        # Perform Wielandt shift
        loss -= 1.0/k_s*prod

        # Begin power iteration
        for i in range(maxits):
            # Check if reach max number of iterations
            if i == maxits - 1:
                raise OpenMCError('Reached maximum iterations in CMFD power '
                                  'iteration solver.')

            # Compute source vector
            s_o = prod.dot(phi_o)

            # Normalize source vector
            s_o /= k_lo

            # Compute new flux with C++ solver
            innerits = openmc.lib._dll.openmc_run_linsolver(loss.data, s_o,
                                                            phi_n, toli)

            # Compute new source vector
            s_n = prod.dot(phi_n)

            # Compute new shifted eigenvalue
            k_ln = np.sum(s_n) / np.sum(s_o)

            # Compute new eigenvalue
            k_n = 1.0/(1.0/k_ln + 1.0/k_s)

            # Renormalize the old source
            s_o *= k_lo

            # Check convergence
            iconv, norm_n = self._check_convergence(s_n, s_o, k_n, k_o, i+1,
                                                    innerits)

            # If converged, calculate dominance ratio and break from loop
            if iconv:
                dom = norm_n / norm_o
                return phi_n, k_n, dom

            # Record old values if not converged
            phi_o = phi_n
            k_o = k_n
            k_lo = k_ln
            norm_o = norm_n

            # Update tolerance for inner iterations
            toli = max(atoli, rtoli*norm_n)

    def _check_convergence(self, s_n, s_o, k_n, k_o, iteration, innerits):
        """Checks the convergence of the CMFD problem

        Parameters
        ----------
        s_n : numpy.ndarray
            Source vector from current iteration
        s_o : numpy.ndarray
            Source vector from previous iteration
        k_n : float
            K-effective  from current iteration
        k_o : float
            K-effective from previous iteration
        iteration : int
            Iteration number
        innerits : int
            Number of iterations required for convergence in inner GS loop

        Returns
        -------
        iconv : bool
            Whether the power iteration has reached convergence
        serr : float
            Error in source from previous iteration to current iteration, used
            for dominance ratio calculations

        """
        # Calculate error in keff
        kerr = abs(k_o - k_n) / k_n

        # Calculate max error in source
        with np.errstate(divide='ignore', invalid='ignore'):
            serr = np.sqrt(np.sum(np.where(s_n > 0, ((s_n-s_o) / s_n)**2, 0))
                           / len(s_n))

        # Check for convergence
        iconv = kerr < self._cmfd_ktol and serr < self._stol

        # Print out to user
        if self._power_monitor and openmc.lib.master():
            print('{:8s}{:20s}{:25s}{:s}{:s}'.format(
                ' {:d}:'.format(iteration),
                'k-eff: {:0.8f}'.format(k_n),
                'k-error:  {:.5e}'.format(kerr),
                'src-error:  {:.5e}'.format(serr),
                '  {:d}'.format(innerits)
            ), flush=True)

        return iconv, serr

    def _set_coremap(self):
        """Sets the core mapping information. All regions marked with zero
        are set to CMFD_NOACCEL, while all regions marked with 1 are set to a
        unique index that maps each fuel region to a row number when building
        CMFD matrices

        """
        # Set number of accelerated regions in problem. This will be related to
        # the dimension of CMFD matrices
        self._mat_dim = np.sum(self._coremap)

        # Define coremap as cumulative sum over accelerated regions,
        # otherwise set value to _CMFD_NOACCEL
        self._coremap = np.where(self._coremap == 0, _CMFD_NOACCEL,
                                 np.cumsum(self._coremap) - 1)

        # Reshape coremap to three dimensional array
        # Indices of coremap in user input switched in x and z axes
        nx, ny, nz = self._indices[:3]
        self._coremap = self._coremap.reshape(nz, ny, nx)
        self._coremap = np.swapaxes(self._coremap, 0, 2)

    def _compute_xs(self):
        """Takes CMFD tallies from OpenMC and computes macroscopic cross
        sections, flux, and diffusion coefficients for each mesh cell using
        a tally window scheme

        """
        # Update window size for expanding window if necessary
        num_cmfd_batches = openmc.lib.current_batch() - self._tally_begin + 1
        if (self._window_type == 'expanding' and
                num_cmfd_batches == self._window_size * 2):
            self._window_size *= 2

        # Discard tallies from oldest batch if window limit reached
        tally_windows = self._flux_rate.shape[-1] + 1
        if tally_windows > self._window_size:
            self._flux_rate = self._flux_rate[...,1:]
            self._total_rate = self._total_rate[...,1:]
            self._p1scatt_rate = self._p1scatt_rate[...,1:]
            self._scatt_rate = self._scatt_rate[...,1:]
            self._nfiss_rate = self._nfiss_rate[...,1:]
            self._current_rate = self._current_rate[...,1:]
            self._openmc_src_rate = self._openmc_src_rate[...,1:]
            tally_windows -= 1

        # Extract spatial and energy indices
        nx, ny, nz, ng = self._indices

        # Get tallies in-memory
        tallies = openmc.lib.tallies

        # Set conditional numpy array as boolean vector based on coremap
        is_accel = self._coremap != _CMFD_NOACCEL

        # Get flux from CMFD tally 0
        tally_id = self._tally_ids[0]
        flux = tallies[tally_id].results[:,0,1]

        # Define target tally reshape dimensions. This defines how openmc
        # tallies are ordered by dimension
        target_tally_shape = [nz, ny, nx, ng, 1]

        # Reshape flux array to target shape. Swap x and z axes so that
        # flux shape is now [nx, ny, nz, ng, 1]
        reshape_flux = np.swapaxes(flux.reshape(target_tally_shape), 0, 2)

        # Flip energy axis as tally results are given in reverse order of
        # energy group
        reshape_flux = np.flip(reshape_flux, axis=3)

        # Bank flux to flux_rate
        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.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 reaction rate (rr) from CMFD tally 0
        totalrr = tallies[tally_id].results[:,1,1]

        # Reshape total reaction rate array to target shape. Swap x and z axes
        # so that shape is now [nx, ny, nz, ng, 1]
        reshape_totalrr = np.swapaxes(totalrr.reshape(target_tally_shape),
                                      0, 2)

        # Total reaction rate is flipped in energy axis as tally results are
        # given in reverse order of energy group
        reshape_totalrr = np.flip(reshape_totalrr, axis=3)

        # Bank total reaction rate to total_rate
        self._total_rate = np.append(self._total_rate, reshape_totalrr,
                                     axis=4)

        # Compute total xs as aggregate of banked total_rate over tally window
        # divided by flux
        self._totalxs = np.divide(np.sum(self._total_rate, axis=4),
                                  self._flux, where=self._flux > 0,
                                  out=np.zeros_like(self._totalxs))

        # Get scattering rr from CMFD tally 1
        # flux is repeated to account for extra dimensionality of scattering xs
        tally_id = self._tally_ids[1]
        scattrr = tallies[tally_id].results[:,0,1]

        # Define target tally reshape dimensions for xs with incoming
        # and outgoing energies
        target_tally_shape = [nz, ny, nx, ng, ng, 1]

        # Reshape scattrr array to target shape. Swap x and z axes so that
        # shape is now [nx, ny, nz, ng, ng, 1]
        reshape_scattrr = np.swapaxes(scattrr.reshape(target_tally_shape),
                                      0, 2)

        # Scattering rr is flipped in both incoming and outgoing energy axes
        # as tally results are given in reverse order of energy group
        reshape_scattrr = np.flip(reshape_scattrr, axis=3)
        reshape_scattrr = np.flip(reshape_scattrr, axis=4)

        # Bank scattering rr to scatt_rate
        self._scatt_rate = np.append(self._scatt_rate, reshape_scattrr,
                                     axis=5)

        # Compute scattering xs as aggregate of banked scatt_rate over tally
        # window divided by flux. Flux dimensionality increased to account for
        # extra dimensionality of scattering xs
        extended_flux = self._flux[:,:,:,:,np.newaxis]
        self._scattxs = np.divide(np.sum(self._scatt_rate, axis=5),
                                  extended_flux, where=extended_flux > 0,
                                  out=np.zeros_like(self._scattxs))

        # Get nu-fission rr from CMFD tally 1
        nfissrr = tallies[tally_id].results[:,1,1]
        num_realizations = tallies[tally_id].num_realizations

        # Reshape nfissrr array to target shape. Swap x and z axes so that
        # shape is now [nx, ny, nz, ng, ng, 1]
        reshape_nfissrr = np.swapaxes(nfissrr.reshape(target_tally_shape),
                                      0, 2)

        # Nu-fission rr is flipped in both incoming and outgoing energy axes
        # as tally results are given in reverse order of energy group
        reshape_nfissrr = np.flip(reshape_nfissrr, axis=3)
        reshape_nfissrr = np.flip(reshape_nfissrr, axis=4)

        # Bank nu-fission rr to nfiss_rate
        self._nfiss_rate = np.append(self._nfiss_rate, reshape_nfissrr,
                                     axis=5)

        # Compute nu-fission xs as aggregate of banked nfiss_rate over tally
        # window divided by flux. Flux dimensionality increased to account for
        # extra dimensionality of nu-fission xs
        self._nfissxs = np.divide(np.sum(self._nfiss_rate, axis=5),
                                  extended_flux, where=extended_flux > 0,
                                  out=np.zeros_like(self._nfissxs))

        # Openmc source distribution is sum of nu-fission rr in incoming
        # energies
        openmc_src = np.sum(reshape_nfissrr, axis=3)

        # Bank OpenMC source distribution from current batch to
        # openmc_src_rate
        self._openmc_src_rate = np.append(self._openmc_src_rate, openmc_src,
                                          axis=4)

        # Compute source distribution over entire tally window
        self._openmc_src = np.sum(self._openmc_src_rate, axis=4)

        # Compute k_eff from source distribution
        self._keff_bal = (np.sum(self._openmc_src) / num_realizations /
                          tally_windows)

        # Normalize openmc source distribution
        self._openmc_src /= np.sum(self._openmc_src) * self._norm

        # Get surface currents from CMFD tally 2
        tally_id = self._tally_ids[2]
        current = tallies[tally_id].results[:,0,1]

        # Define target tally reshape dimensions for current
        target_tally_shape = [nz, ny, nx, 12, ng, 1]

        # Reshape current array to target shape. Swap x and z axes so that
        # shape is now [nx, ny, nz, 12, ng, 1]
        reshape_current = np.swapaxes(current.reshape(target_tally_shape),
                                      0, 2)

        # Current is flipped in energy axis as tally results are given in
        # reverse order of energy group
        reshape_current = np.flip(reshape_current, axis=4)

        # Bank current to current_rate
        self._current_rate = np.append(self._current_rate, reshape_current,
                                       axis=5)

        # Compute current as aggregate of banked current_rate over tally window
        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]
        p1scattrr = tallies[tally_id].results[:,0,1]

        # Define target tally reshape dimensions for p1 scatter tally
        target_tally_shape = [nz, ny, nx, 2, ng, 1]

        # Reshape and extract only p1 data from tally results as there is
        # no need for p0 data
        reshape_p1scattrr = np.swapaxes(p1scattrr.reshape(target_tally_shape),
                                        0, 2)[:,:,:,1,:,:]

        # p1-scatter rr is flipped in energy axis as tally results are given in
        # reverse order of energy group
        reshape_p1scattrr = np.flip(reshape_p1scattrr, axis=3)

        # Bank p1-scatter rr to p1scatt_rate
        self._p1scatt_rate = np.append(self._p1scatt_rate, reshape_p1scattrr,
                                       axis=4)

        # Compute p1-scatter xs as aggregate of banked p1scatt_rate over tally
        # window divided by flux
        self._p1scattxs = np.divide(np.sum(self._p1scatt_rate, axis=4),
                                    self._flux, where=self._flux > 0,
                                    out=np.zeros_like(self._p1scattxs))

        if self._set_reference_params:
            # Set diffusion coefficients based on reference value
            self._diffcof = np.where(self._flux > 0,
                                     self._ref_d[None, None, None, :], 0.0)
        else:
            # Calculate and store diffusion coefficient
            with np.errstate(divide='ignore', invalid='ignore'):
                self._diffcof = np.where(self._flux > 0, 1.0 / (3.0 *
                                         (self._totalxs-self._p1scattxs)), 0.)

    def _compute_effective_downscatter(self):
        """Changes downscatter rate for zero upscatter"""
        # Extract energy index
        ng = self._indices[3]

        # Return if not two groups
        if ng != 2:
            return

        # Extract cross sections and flux for each group
        flux1 = self._flux[:,:,:,0]
        flux2 = self._flux[:,:,:,1]
        sigt1 = self._totalxs[:,:,:,0]
        sigt2 = self._totalxs[:,:,:,1]

        # First energy index is incoming energy, second is outgoing energy
        sigs11 = self._scattxs[:,:,:,0,0]
        sigs21 = self._scattxs[:,:,:,1,0]
        sigs12 = self._scattxs[:,:,:,0,1]
        sigs22 = self._scattxs[:,:,:,1,1]

        # Compute absorption xs
        siga1 = sigt1 - sigs11 - sigs12
        siga2 = sigt2 - sigs22 - sigs21

        # Compute effective downscatter XS
        sigs12_eff = sigs12 - sigs21 * np.divide(flux2, flux1,
                                                 where=flux1 > 0,
                                                 out=np.zeros_like(flux2))

        # Recompute total cross sections and record
        self._totalxs[:,:,:,0] = siga1 + sigs11 + sigs12_eff
        self._totalxs[:,:,:,1] = siga2 + sigs22

        # Record effective dowmscatter xs
        self._scattxs[:,:,:,0,1] = sigs12_eff

        # Zero out upscatter cross section
        self._scattxs[:,:,:,1,0] = 0.0

    def _neutron_balance(self):
        """Computes the RMS neutron balance over the CMFD mesh"""
        # Extract energy indices
        ng = self._indices[3]

        # Get number of accelerated regions
        num_accel = self._mat_dim

        # Get openmc k-effective
        keff = openmc.lib.keff()[0]

        # Define leakage in each mesh cell and energy group
        leakage = (((self._current[:,:,:,_CURRENTS['out_right'],:] -
                   self._current[:,:,:,_CURRENTS['in_right'],:]) -
                   (self._current[:,:,:,_CURRENTS['in_left'],:] -
                   self._current[:,:,:,_CURRENTS['out_left'],:])) +
                   ((self._current[:,:,:,_CURRENTS['out_front'],:] -
                    self._current[:,:,:,_CURRENTS['in_front'],:]) -
                   (self._current[:,:,:,_CURRENTS['in_back'],:] -
                    self._current[:,:,:,_CURRENTS['out_back'],:])) +
                   ((self._current[:,:,:,_CURRENTS['out_top'],:] -
                    self._current[:,:,:,_CURRENTS['in_top'],:]) -
                   (self._current[:,:,:,_CURRENTS['in_bottom'],:] -
                    self._current[:,:,:,_CURRENTS['out_bottom'],:])))

        # Compute total rr
        interactions = self._totalxs * self._flux

        # Compute scattering rr by broadcasting flux in outgoing energy and
        # summing over incoming energy
        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],
                         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,
                                out=np.zeros_like(self._flux))

        # Calculate RMS and record for this batch
        self._balance.append(np.sqrt(
            np.sum(np.multiply(self._resnb, self._resnb)) /
            (ng * num_accel)))

    def _precompute_array_indices(self):
        """Initializes cross section arrays and computes the indices
        used to populate dtilde and dhat

        """
        # Extract spatial 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.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
        is_accel = self._coremap != _CMFD_NOACCEL
        # Logical for determining whether a zero flux "albedo" b.c. should be
        # applied
        x_inds, y_inds, z_inds = np.indices((nx, ny, nz))

        # Define slice equivalent to is_accel[0,:,:]
        slice_x = x_inds[:1,:,:]
        slice_y = y_inds[:1,:,:]
        slice_z = z_inds[:1,:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._first_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                               slice_z[bndry_accel])

        # Define slice equivalent to is_accel[-1,:,:]
        slice_x = x_inds[-1:,:,:]
        slice_y = y_inds[-1:,:,:]
        slice_z = z_inds[-1:,:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._last_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                              slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,0,:]
        slice_x = x_inds[:,:1,:]
        slice_y = y_inds[:,:1,:]
        slice_z = z_inds[:,:1,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._first_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                               slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,-1,:]
        slice_x = x_inds[:,-1:,:]
        slice_y = y_inds[:,-1:,:]
        slice_z = z_inds[:,-1:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._last_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                              slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,:,0]
        slice_x = x_inds[:,:,:1]
        slice_y = y_inds[:,:,:1]
        slice_z = z_inds[:,:,:1]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._first_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                               slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,:,-1]
        slice_x = x_inds[:,:,-1:]
        slice_y = y_inds[:,:,-1:]
        slice_z = z_inds[:,:,-1:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._last_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                              slice_z[bndry_accel])

        # Define slice equivalent to is_accel[1:,:,:]
        slice_x = x_inds[1:,:,:]
        slice_y = y_inds[1:,:,:]
        slice_z = z_inds[1:,:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notfirst_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                  slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:-1,:,:]
        slice_x = x_inds[:-1,:,:]
        slice_y = y_inds[:-1,:,:]
        slice_z = z_inds[:-1,:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notlast_x_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                 slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,1:,:]
        slice_x = x_inds[:,1:,:]
        slice_y = y_inds[:,1:,:]
        slice_z = z_inds[:,1:,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notfirst_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                  slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,:-1,:]
        slice_x = x_inds[:,:-1,:]
        slice_y = y_inds[:,:-1,:]
        slice_z = z_inds[:,:-1,:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notlast_y_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                 slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,:,1:]
        slice_x = x_inds[:,:,1:]
        slice_y = y_inds[:,:,1:]
        slice_z = z_inds[:,:,1:]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notfirst_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                  slice_z[bndry_accel])

        # Define slice equivalent to is_accel[:,:,:-1]
        slice_x = x_inds[:,:,:-1]
        slice_y = y_inds[:,:,:-1]
        slice_z = z_inds[:,:,:-1]
        bndry_accel = is_accel[(slice_x, slice_y, slice_z)]
        self._notlast_z_accel = (slice_x[bndry_accel], slice_y[bndry_accel],
                                 slice_z[bndry_accel])

        # Store logical for whether neighboring cell is reflector region
        # in all directions
        adj_reflector_left = np.roll(self._coremap, 1, axis=0) == _CMFD_NOACCEL
        self._is_adj_ref_left = adj_reflector_left[
                self._notfirst_x_accel + (np.newaxis,)]

        adj_reflector_right = np.roll(self._coremap, -1, axis=0) == \
            _CMFD_NOACCEL
        self._is_adj_ref_right = adj_reflector_right[
                self._notlast_x_accel + (np.newaxis,)]

        adj_reflector_back = np.roll(self._coremap, 1, axis=1) == \
            _CMFD_NOACCEL
        self._is_adj_ref_back = adj_reflector_back[
                self._notfirst_y_accel + (np.newaxis,)]

        adj_reflector_front = np.roll(self._coremap, -1, axis=1) == \
            _CMFD_NOACCEL
        self._is_adj_ref_front = adj_reflector_front[
                self._notlast_y_accel + (np.newaxis,)]

        adj_reflector_bottom = np.roll(self._coremap, 1, axis=2) == \
            _CMFD_NOACCEL
        self._is_adj_ref_bottom = adj_reflector_bottom[
                self._notfirst_z_accel + (np.newaxis,)]

        adj_reflector_top = np.roll(self._coremap, -1, axis=2) == \
            _CMFD_NOACCEL
        self._is_adj_ref_top = adj_reflector_top[
                self._notlast_z_accel + (np.newaxis,)]

    def _precompute_matrix_indices(self):
        """Computes the indices and row/column data used to populate CMFD CSR
        matrices. These indices are used in _build_loss_matrix and
        _build_prod_matrix.

        """
        # Extract energy group indices
        ng = self._indices[3]

        # Shift coremap in all directions to determine whether leakage term
        # should be defined for particular cell in matrix
        coremap_shift_left = np.pad(self._coremap, ((1,0),(0,0),(0,0)),
                                    mode='constant',
                                    constant_values=_CMFD_NOACCEL)[:-1,:,:]

        coremap_shift_right = np.pad(self._coremap, ((0,1),(0,0),(0,0)),
                                     mode='constant',
                                     constant_values=_CMFD_NOACCEL)[1:,:,:]

        coremap_shift_back = np.pad(self._coremap, ((0,0),(1,0),(0,0)),
                                    mode='constant',
                                    constant_values=_CMFD_NOACCEL)[:,:-1,:]

        coremap_shift_front = np.pad(self._coremap, ((0,0),(0,1),(0,0)),
                                     mode='constant',
                                     constant_values=_CMFD_NOACCEL)[:,1:,:]

        coremap_shift_bottom = np.pad(self._coremap, ((0,0),(0,0),(1,0)),
                                      mode='constant',
                                      constant_values=_CMFD_NOACCEL)[:,:,:-1]

        coremap_shift_top = np.pad(self._coremap, ((0,0),(0,0),(0,1)),
                                   mode='constant',
                                   constant_values=_CMFD_NOACCEL)[:,:,1:]

        # Create empty row and column vectors to store for loss matrix
        row = np.array([])
        col = np.array([])

        # Store all indices used to populate production and loss matrix
        is_accel = self._coremap != _CMFD_NOACCEL
        self._accel_idxs = np.where(is_accel)
        self._accel_neig_left_idxs = (np.where(is_accel &
                                      (coremap_shift_left != _CMFD_NOACCEL)))
        self._accel_neig_right_idxs = (np.where(is_accel &
                                       (coremap_shift_right != _CMFD_NOACCEL)))
        self._accel_neig_back_idxs = (np.where(is_accel &
                                      (coremap_shift_back != _CMFD_NOACCEL)))
        self._accel_neig_front_idxs = (np.where(is_accel &
                                       (coremap_shift_front != _CMFD_NOACCEL)))
        self._accel_neig_bot_idxs = (np.where(is_accel &
                                     (coremap_shift_bottom != _CMFD_NOACCEL)))
        self._accel_neig_top_idxs = (np.where(is_accel &
                                     (coremap_shift_top != _CMFD_NOACCEL)))

        for g in range(ng):
            # Extract row and column data of regions where a cell and its
            # neighbor to the left are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_left_idxs]) + g
            idx_y = ng * (coremap_shift_left[self._accel_neig_left_idxs]) + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract row and column data of regions where a cell and its
            # neighbor to the right are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_right_idxs]) + g
            idx_y = ng * (coremap_shift_right[self._accel_neig_right_idxs]) + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract row and column data of regions where a cell and its
            # neighbor to the back are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_back_idxs]) + g
            idx_y = ng * (coremap_shift_back[self._accel_neig_back_idxs]) + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract row and column data of regions where a cell and its
            # neighbor to the front are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_front_idxs]) + g
            idx_y = ng * (coremap_shift_front[self._accel_neig_front_idxs]) + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract row and column data of regions where a cell and its
            # neighbor to the bottom are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_bot_idxs]) + g
            idx_y = ng * (coremap_shift_bottom[self._accel_neig_bot_idxs]) \
                + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract row and column data of regions where a cell and its
            # neighbor to the top are both fuel regions
            idx_x = ng * (self._coremap[self._accel_neig_top_idxs]) + g
            idx_y = ng * (coremap_shift_top[self._accel_neig_top_idxs]) + g
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            # Extract all regions where a cell is a fuel region
            idx_x = ng * (self._coremap[self._accel_idxs]) + g
            idx_y = idx_x
            row = np.append(row, idx_x)
            col = np.append(col, idx_y)

            for h in range(ng):
                if h != g:
                    # Extract all regions where a cell is a fuel region
                    idx_x = ng * (self._coremap[self._accel_idxs]) + g
                    idx_y = ng * (self._coremap[self._accel_idxs]) + h
                    row = np.append(row, idx_x)
                    col = np.append(col, idx_y)

        # Store row and col as rows and columns of production matrix
        self._loss_row = row
        self._loss_col = col

        # Create empty row and column vectors to store for production matrix
        row = np.array([], dtype=int)
        col = np.array([], dtype=int)

        for g in range(ng):
            for h in range(ng):
                # Extract all regions where a cell is a fuel region
                idx_x = ng * (self._coremap[self._accel_idxs]) + g
                idx_y = ng * (self._coremap[self._accel_idxs]) + h
                # Store rows, cols, and data to add to CSR matrix
                row = np.append(row, idx_x)
                col = np.append(col, idx_y)

        # Store row and col as rows and columns of production matrix
        self._prod_row = row
        self._prod_col = col

    def _compute_dtilde(self):
        """Computes the diffusion coupling coefficient using a vectorized numpy
        approach. Aggregate values for the dtilde multidimensional array are
        populated by first defining values on the problem boundary, and then
        for all other regions. For indices not lying on a boundary, dtilde
        values are distinguished between regions that neighbor a reflector
        region and regions that don't neighbor a reflector

        """
        # Logical for determining whether a zero flux "albedo" b.c. should be
        # applied
        is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT

        # Define dtilde at left surface for all mesh cells on left boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._first_x_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dx = self._hxyz[boundary + (np.newaxis, 0)]
        if is_zero_flux_alb[0]:
            self._dtilde[boundary_grps + (0,)] = 2.0 * D / dx
        else:
            alb = self._albedo[0]
            self._dtilde[boundary_grps + (0,)] = ((2.0 * D * (1.0 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dx))

        # Define dtilde at right surface for all mesh cells on right boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._last_x_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dx = self._hxyz[boundary + (np.newaxis, 0)]
        if is_zero_flux_alb[1]:
            self._dtilde[boundary_grps + (1,)] = 2.0 * D / dx
        else:
            alb = self._albedo[1]
            self._dtilde[boundary_grps + (1,)] = ((2.0 * D * (1.0 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dx))

        # Define dtilde at back surface for all mesh cells on back boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._first_y_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dy = self._hxyz[boundary + (np.newaxis, 1)]
        if is_zero_flux_alb[2]:
            self._dtilde[boundary_grps + (2,)] = 2.0 * D / dy
        else:
            alb = self._albedo[2]
            self._dtilde[boundary_grps + (2,)] = ((2.0 * D * (1.0 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dy))

        # Define dtilde at front surface for all mesh cells on front boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._last_y_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dy = self._hxyz[boundary + (np.newaxis, 1)]
        if is_zero_flux_alb[3]:
            self._dtilde[boundary_grps + (3,)] = 2.0 * D / dy
        else:
            alb = self._albedo[3]
            self._dtilde[boundary_grps + (3,)] = ((2.0 * D * (1.0 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dy))

        # Define dtilde at bottom surface for all mesh cells on bottom boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._first_z_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dz = self._hxyz[boundary + (np.newaxis, 2)]
        if is_zero_flux_alb[4]:
            self._dtilde[boundary_grps + (4,)] = 2.0 * D / dz
        else:
            alb = self._albedo[4]
            self._dtilde[boundary_grps + (4,)] = ((2.0 * D * (1.0 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dz))

        # Define dtilde at top surface for all mesh cells on top boundary
        # Separate between zero flux b.c. and alebdo b.c.
        boundary = self._last_z_accel
        boundary_grps = boundary + (slice(None),)

        D = self._diffcof[boundary_grps]
        dz = self._hxyz[boundary + (np.newaxis, 2)]
        if is_zero_flux_alb[5]:
            self._dtilde[boundary_grps + (5,)] = 2.0 * D / dz
        else:
            alb = self._albedo[5]
            self._dtilde[boundary_grps + (5,)] = ((2.0 * D * (1 - alb))
                                                  / (4.0 * D * (1.0 + alb) +
                                                  (1.0 - alb) * dz))

        # Define reflector albedo for all cells on the left surface, in case
        # a cell borders a reflector region on the left
        current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:]
        current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:]
        ref_albedo = np.divide(current_in_left, current_out_left,
                               where=current_out_left > 1.0e-10,
                               out=np.ones_like(current_out_left))

        # Diffusion coefficient of neighbor to left
        neig_dc = np.roll(self._diffcof, 1, axis=0)
        # Cell dimensions of neighbor to left
        neig_hxyz = np.roll(self._hxyz, 1, axis=0)

        # Define dtilde at left surface for all mesh cells not on left boundary
        # Dtilde is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notfirst_x_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dx = self._hxyz[boundary + (np.newaxis, 0)]
        neig_D = neig_dc[boundary_grps]
        neig_dx = neig_hxyz[boundary + (np.newaxis, 0)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_left
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx),
                          (2.0 * D * neig_D) / (neig_dx * D + dx * neig_D))
        self._dtilde[boundary_grps + (0,)] = dtilde

        # Define reflector albedo for all cells on the right surface, in case
        # a cell borders a reflector region on the right
        current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:]
        current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:]
        ref_albedo = np.divide(current_in_right, current_out_right,
                               where=current_out_right > 1.0e-10,
                               out=np.ones_like(current_out_right))

        # Diffusion coefficient of neighbor to right
        neig_dc = np.roll(self._diffcof, -1, axis=0)
        # Cell dimensions of neighbor to right
        neig_hxyz = np.roll(self._hxyz, -1, axis=0)

        # Define dtilde at right surface for all mesh cells not on right
        # boundary. Dtilde is defined differently for regions that do and don't
        # neighbor reflector regions
        boundary = self._notlast_x_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dx = self._hxyz[boundary + (np.newaxis, 0)]
        neig_D = neig_dc[boundary_grps]
        neig_dx = neig_hxyz[boundary + (np.newaxis, 0)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_right
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dx),
                          (2.0 * D * neig_D) / (neig_dx * D + dx * neig_D))
        self._dtilde[boundary_grps + (1,)] = dtilde

        # Define reflector albedo for all cells on the back surface, in case
        # a cell borders a reflector region on the back
        current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:]
        current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:]
        ref_albedo = np.divide(current_in_back, current_out_back,
                               where=current_out_back > 1.0e-10,
                               out=np.ones_like(current_out_back))

        # Diffusion coefficient of neighbor to back
        neig_dc = np.roll(self._diffcof, 1, axis=1)
        # Cell dimensions of neighbor to back
        neig_hxyz = np.roll(self._hxyz, 1, axis=1)

        # Define dtilde at back surface for all mesh cells not on back boundary
        # Dtilde is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notfirst_y_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dy = self._hxyz[boundary + (np.newaxis, 1)]
        neig_D = neig_dc[boundary_grps]
        neig_dy = neig_hxyz[boundary + (np.newaxis, 1)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_back
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy),
                          (2.0 * D * neig_D) / (neig_dy * D + dy * neig_D))
        self._dtilde[boundary_grps + (2,)] = dtilde

        # Define reflector albedo for all cells on the front surface, in case
        # a cell borders a reflector region in the front
        current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:]
        current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:]
        ref_albedo = np.divide(current_in_front, current_out_front,
                               where=current_out_front > 1.0e-10,
                               out=np.ones_like(current_out_front))

        # Diffusion coefficient of neighbor to front
        neig_dc = np.roll(self._diffcof, -1, axis=1)
        # Cell dimensions of neighbor to front
        neig_hxyz = np.roll(self._hxyz, -1, axis=1)

        # Define dtilde at front surface for all mesh cells not on front
        # boundary. Dtilde is defined differently for regions that do and don't
        # neighbor reflector regions
        boundary = self._notlast_y_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dy = self._hxyz[boundary + (np.newaxis, 1)]
        neig_D = neig_dc[boundary_grps]
        neig_dy = neig_hxyz[boundary + (np.newaxis, 1)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_front
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dy),
                          (2.0 * D * neig_D) / (neig_dy * D + dy * neig_D))
        self._dtilde[boundary_grps + (3,)] = dtilde

        # Define reflector albedo for all cells on the bottom surface, in case
        # a cell borders a reflector region on the bottom
        current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:]
        current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:]
        ref_albedo = np.divide(current_in_bottom, current_out_bottom,
                               where=current_out_bottom > 1.0e-10,
                               out=np.ones_like(current_out_bottom))

        # Diffusion coefficient of neighbor to bottom
        neig_dc = np.roll(self._diffcof, 1, axis=2)
        # Cell dimensions of neighbor to bottom
        neig_hxyz = np.roll(self._hxyz, 1, axis=2)

        # Define dtilde at bottom surface for all mesh cells not on bottom
        # boundary. Dtilde is defined differently for regions that do and don't
        # neighbor reflector regions
        boundary = self._notfirst_z_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dz = self._hxyz[boundary + (np.newaxis, 2)]
        neig_D = neig_dc[boundary_grps]
        neig_dz = neig_hxyz[boundary + (np.newaxis, 2)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_bottom
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz),
                          (2.0 * D * neig_D) / (neig_dz * D + dz * neig_D))
        self._dtilde[boundary_grps + (4,)] = dtilde

        # Define reflector albedo for all cells on the top surface, in case
        # a cell borders a reflector region on the top
        current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:]
        current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:]
        ref_albedo = np.divide(current_in_top, current_out_top,
                               where=current_out_top > 1.0e-10,
                               out=np.ones_like(current_out_top))

        # Diffusion coefficient of neighbor to top
        neig_dc = np.roll(self._diffcof, -1, axis=2)
        # Cell dimensions of neighbor to top
        neig_hxyz = np.roll(self._hxyz, -1, axis=2)

        # Define dtilde at top surface for all mesh cells not on top boundary
        # Dtilde is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notlast_z_accel
        boundary_grps = boundary + (slice(None),)
        D = self._diffcof[boundary_grps]
        dz = self._hxyz[boundary + (np.newaxis, 2)]
        neig_D = neig_dc[boundary_grps]
        neig_dz = neig_hxyz[boundary + (np.newaxis, 2)]
        alb = ref_albedo[boundary_grps]
        is_adj_ref = self._is_adj_ref_top
        dtilde = np.where(is_adj_ref, (2.0 * D * (1.0 - alb)) /
                          (4.0 * D * (1.0 + alb) + (1.0 - alb) * dz),
                          (2.0 * D * neig_D) / (neig_dz * D + dz * neig_D))
        self._dtilde[boundary_grps + (5,)] = dtilde

    def _compute_dhat(self):
        """Computes the nonlinear coupling coefficient using a vectorized numpy
        approach. Aggregate values for the dhat multidimensional array are
        populated by first defining values on the problem boundary, and then
        for all other regions. For indices not lying by a boundary, dhat values
        are distinguished between regions that neighbor a reflector region and
        regions that don't neighbor a reflector

        """
        # Define current in each direction
        current_in_left = self._current[:,:,:,_CURRENTS['in_left'],:]
        current_out_left = self._current[:,:,:,_CURRENTS['out_left'],:]
        current_in_right = self._current[:,:,:,_CURRENTS['in_right'],:]
        current_out_right = self._current[:,:,:,_CURRENTS['out_right'],:]
        current_in_back = self._current[:,:,:,_CURRENTS['in_back'],:]
        current_out_back = self._current[:,:,:,_CURRENTS['out_back'],:]
        current_in_front = self._current[:,:,:,_CURRENTS['in_front'],:]
        current_out_front = self._current[:,:,:,_CURRENTS['out_front'],:]
        current_in_bottom = self._current[:,:,:,_CURRENTS['in_bottom'],:]
        current_out_bottom = self._current[:,:,:,_CURRENTS['out_bottom'],:]
        current_in_top = self._current[:,:,:,_CURRENTS['in_top'],:]
        current_out_top = self._current[:,:,:,_CURRENTS['out_top'],:]

        dx = self._hxyz[:,:,:,np.newaxis,0]
        dy = self._hxyz[:,:,:,np.newaxis,1]
        dz = self._hxyz[:,:,:,np.newaxis,2]
        dxdydz = np.prod(self._hxyz, axis=3)[:,:,:,np.newaxis]

        # Define net current on each face
        net_current_left = (current_in_left - current_out_left) / dxdydz * dx
        net_current_right = (current_out_right - current_in_right) / dxdydz * \
            dx
        net_current_back = (current_in_back - current_out_back) / dxdydz * dy
        net_current_front = (current_out_front - current_in_front) / dxdydz * \
            dy
        net_current_bottom = (current_in_bottom - current_out_bottom) / \
            dxdydz * dz
        net_current_top = (current_out_top - current_in_top) / dxdydz * dz

        # Define flux in each cell
        cell_flux = self._flux / dxdydz

        # Define dhat at left surface for all mesh cells on left boundary
        boundary = self._first_x_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_left[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 0)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (0,)] = (net_current + dtilde * flux) / flux

        # Define dhat at right surface for all mesh cells on right boundary
        boundary = self._last_x_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_right[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 1)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (1,)] = (net_current - dtilde * flux) / flux

        # Define dhat at back surface for all mesh cells on back boundary
        boundary = self._first_y_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_back[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 2)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (2,)] = (net_current + dtilde * flux) / flux

        # Define dhat at front surface for all mesh cells on front boundary
        boundary = self._last_y_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_front[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 3)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (3,)] = (net_current - dtilde * flux) / flux

        # Define dhat at bottom surface for all mesh cells on bottom boundary
        boundary = self._first_z_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_bottom[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 4)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (4,)] = (net_current + dtilde * flux) / flux

        # Define dhat at top surface for all mesh cells on top boundary
        boundary = self._last_z_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_top[boundary_grps]
        dtilde = self._dtilde[boundary + (slice(None), 5)]
        flux = cell_flux[boundary_grps]
        self._dhat[boundary_grps + (5,)] = (net_current - dtilde * flux) / flux

        # Cell flux of neighbor to left
        neig_flux = np.roll(self._flux, 1, axis=0) / dxdydz

        # Define dhat at left surface for all mesh cells not on left boundary
        # Dhat is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notfirst_x_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_left[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (0,)]
        flux = cell_flux[boundary_grps]
        flux_left = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_left
        dhat = np.where(is_adj_ref, (net_current + dtilde * flux) / flux,
                        (net_current - dtilde * (flux_left - flux)) /
                        (flux_left + flux))
        self._dhat[boundary_grps + (0,)] = dhat

        # Cell flux of neighbor to right
        neig_flux = np.roll(self._flux, -1, axis=0) / dxdydz

        # Define dhat at right surface for all mesh cells not on right boundary
        # Dhat is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notlast_x_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_right[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (1,)]
        flux = cell_flux[boundary_grps]
        flux_right = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_right
        dhat = np.where(is_adj_ref, (net_current - dtilde * flux) / flux,
                        (net_current + dtilde * (flux_right - flux)) /
                        (flux_right + flux))
        self._dhat[boundary_grps + (1,)] = dhat

        # Cell flux of neighbor to back
        neig_flux = np.roll(self._flux, 1, axis=1) / dxdydz

        # Define dhat at back surface for all mesh cells not on back boundary
        # Dhat is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notfirst_y_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_back[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (2,)]
        flux = cell_flux[boundary_grps]
        flux_back = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_back
        dhat = np.where(is_adj_ref, (net_current + dtilde * flux) / flux,
                        (net_current - dtilde * (flux_back - flux)) /
                        (flux_back + flux))
        self._dhat[boundary_grps + (2,)] = dhat

        # Cell flux of neighbor to front
        neig_flux = np.roll(self._flux, -1, axis=1) / dxdydz

        # Define dhat at front surface for all mesh cells not on front boundary
        # Dhat is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notlast_y_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_front[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (3,)]
        flux = cell_flux[boundary_grps]
        flux_front = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_front
        dhat = np.where(is_adj_ref, (net_current - dtilde * flux) / flux,
                        (net_current + dtilde * (flux_front - flux)) /
                        (flux_front + flux))
        self._dhat[boundary_grps + (3,)] = dhat

        # Cell flux of neighbor to bottom
        neig_flux = np.roll(self._flux, 1, axis=2) / dxdydz

        # Define dhat at bottom surface for all mesh cells not on bottom
        # boundary. Dhat is defined differently for regions that do and don't
        # neighbor reflector regions
        boundary = self._notfirst_z_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_bottom[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (4,)]
        flux = cell_flux[boundary_grps]
        flux_bottom = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_bottom
        dhat = np.where(is_adj_ref, (net_current + dtilde * flux) / flux,
                        (net_current - dtilde * (flux_bottom - flux)) /
                        (flux_bottom + flux))
        self._dhat[boundary_grps + (4,)] = dhat

        # Cell flux of neighbor to top
        neig_flux = np.roll(self._flux, -1, axis=2) / dxdydz

        # Define dhat at top surface for all mesh cells not on top boundary
        # Dhat is defined differently for regions that do and don't neighbor
        # reflector regions
        boundary = self._notlast_z_accel
        boundary_grps = boundary + (slice(None),)
        net_current = net_current_top[boundary_grps]
        dtilde = self._dtilde[boundary_grps + (5,)]
        flux = cell_flux[boundary_grps]
        flux_top = neig_flux[boundary_grps]
        is_adj_ref = self._is_adj_ref_top
        dhat = np.where(is_adj_ref, (net_current - dtilde * flux) / flux,
                        (net_current + dtilde * (flux_top - flux)) /
                        (flux_top + flux))
        self._dhat[boundary_grps + (5,)] = dhat

    def _create_cmfd_tally(self):
        """Creates all tallies in-memory that are used to solve CMFD problem"""
        # Create Mesh object based on CMFDMesh mesh_type, stored internally
        if self._mesh.mesh_type == 'regular':
            cmfd_mesh = openmc.lib.RegularMesh()
            # Set dimension and parameters of mesh object
            cmfd_mesh.dimension = self._mesh.dimension
            cmfd_mesh.set_parameters(lower_left=self._mesh.lower_left,
                                     upper_right=self._mesh.upper_right,
                                     width=self._mesh.width)
        elif self._mesh.mesh_type == 'rectilinear':
            cmfd_mesh = openmc.lib.RectilinearMesh()
            # Set grid of mesh object
            x_grid, y_grid, z_grid = self._mesh.grid
            cmfd_mesh.set_grid(x_grid, y_grid, z_grid)

        # Store id of mesh object
        self._mesh_id = cmfd_mesh.id

        # Create mesh Filter object, stored internally
        mesh_filter = openmc.lib.MeshFilter()
        # Set mesh for Mesh Filter
        mesh_filter.mesh = cmfd_mesh

        # Set up energy filters, if applicable
        if self._energy_filters:
            # Create Energy Filter object, stored internally
            energy_filter = openmc.lib.EnergyFilter()
            # Set bins for Energy Filter
            energy_filter.bins = self._egrid

            # Create Energy Out Filter object, stored internally
            energyout_filter = openmc.lib.EnergyoutFilter()
            # Set bins for Energy Filter
            energyout_filter.bins = self._egrid

        # Create Mesh Surface Filter object, stored internally
        meshsurface_filter = openmc.lib.MeshSurfaceFilter()
        # Set mesh for Mesh Surface Filter
        meshsurface_filter.mesh = cmfd_mesh

        # Create Legendre Filter object, stored internally
        legendre_filter = openmc.lib.LegendreFilter()
        # Set order for Legendre Filter
        legendre_filter.order = 1

        # Create CMFD tallies, stored internally
        n_tallies = 4
        self._tally_ids = []
        for i in range(n_tallies):
            cmfd_tally = openmc.lib.Tally()
            # Set nuclide bins
            cmfd_tally.nuclides = ['total']
            self._tally_ids.append(cmfd_tally.id)

            # Set attributes of CMFD flux, total tally
            if i == 0:
                # Set filters for tally
                if self._energy_filters:
                    cmfd_tally.filters = [mesh_filter, energy_filter]
                else:
                    cmfd_tally.filters = [mesh_filter]
                # Set scores, type, and estimator for tally
                cmfd_tally.scores = ['flux', 'total']
                cmfd_tally.type = 'volume'
                cmfd_tally.estimator = 'analog'

            # Set attributes of CMFD neutron production tally
            elif i == 1:
                # Set filters for tally
                if self._energy_filters:
                    cmfd_tally.filters = [mesh_filter, energy_filter,
                                          energyout_filter]
                else:
                    cmfd_tally.filters = [mesh_filter]
                # Set scores, type, and estimator for tally
                cmfd_tally.scores = ['nu-scatter', 'nu-fission']
                cmfd_tally.type = 'volume'
                cmfd_tally.estimator = 'analog'

            # Set attributes of CMFD surface current tally
            elif i == 2:
                # Set filters for tally
                if self._energy_filters:
                    cmfd_tally.filters = [meshsurface_filter, energy_filter]
                else:
                    cmfd_tally.filters = [meshsurface_filter]
                # Set scores, type, and estimator for tally
                cmfd_tally.scores = ['current']
                cmfd_tally.type = 'mesh-surface'
                cmfd_tally.estimator = 'analog'

            # Set attributes of CMFD P1 scatter tally
            elif i == 3:
                # Set filters for tally
                if self._energy_filters:
                    cmfd_tally.filters = [mesh_filter, legendre_filter,
                                          energy_filter]
                else:
                    cmfd_tally.filters = [mesh_filter, legendre_filter]
                # Set scores for tally
                cmfd_tally.scores = ['scatter']
                cmfd_tally.type = 'volume'
                cmfd_tally.estimator = 'analog'

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