crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1124 from shikhar413/cmfd-capi

Browse source 
Move CMFD to Python through C API - Replace Fortran, update docs and test suite
This commit is contained in:
Sterling Harper 2018-12-23 21:31:10 -05:00 • committed by GitHub
commit 4517896bdb
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
60 changed files with 4491 additions and 4696 deletions
Download patch file Download diff file Expand all files Collapse all files

9
CMakeLists.txt
View file

@ -309,13 +309,6 @@ add_library(libopenmc SHARED
src/bank_header.F90
src/api.F90
src/dagmc_header.F90
src/cmfd_data.F90
src/cmfd_execute.F90
src/cmfd_header.F90
src/cmfd_input.F90
src/cmfd_loss_operator.F90
src/cmfd_prod_operator.F90
src/cmfd_solver.F90
src/constants.F90
src/dict_header.F90
src/eigenvalue.F90
@ -383,7 +376,7 @@ add_library(libopenmc SHARED
src/bank.cpp
src/dagmc.cpp
src/cell.cpp
src/cmfd_execute.cpp
src/cmfd_solver.cpp
src/cross_sections.cpp
src/distribution.cpp
src/distribution_angle.cpp

221
docs/source/io_formats/cmfd.rst
View file

@ -1,221 +0,0 @@
.. _io_cmfd:
==============================
CMFD Specification -- cmfd.xml
==============================
Coarse mesh finite difference acceleration method has been implemented in
OpenMC. Currently, it allows users to accelerate fission source convergence
during inactive neutron batches. To run CMFD, the ``<run_cmfd>`` element in
``settings.xml`` should be set to "true".
-------------------
``<begin>`` Element
-------------------
The ``<begin>`` element controls what batch CMFD calculations should begin.
*Default*: 1
------------------------
``<dhat_reset>`` Element
------------------------
The ``<dhat_reset>`` element controls whether :math:`\widehat{D}` nonlinear
CMFD parameters should be reset to zero before solving CMFD eigenproblem.
It can be turned on with "true" and off with "false".
*Default*: false
---------------------
``<display>`` Element
---------------------
The ``<display>`` element sets one additional CMFD output column. Options are:
* "balance" - prints the RMS [%] of the resdiual from the neutron balance
equation on CMFD tallies.
* "dominance" - prints the estimated dominance ratio from the CMFD iterations.
**This will only work for power iteration eigensolver**.
* "entropy" - prints the *entropy* of the CMFD predicted fission source.
**Can only be used if OpenMC entropy is active as well**.
* "source" - prints the RMS [%] between the OpenMC fission source and CMFD
fission source.
*Default*: balance
-------------------------
``<downscatter>`` Element
-------------------------
The ``<downscatter>`` element controls whether an effective downscatter cross
section should be used when using 2-group CMFD. It can be turned on with "true"
and off with "false".
*Default*: false
----------------------
``<feedback>`` Element
----------------------
The ``<feedback>`` element controls whether or not the CMFD diffusion result is
used to adjust the weight of fission source neutrons on the next OpenMC batch.
It can be turned on with "true" and off with "false".
*Default*: false
------------------------------------
``<gauss_seidel_tolerance>`` Element
------------------------------------
The ``<gauss_seidel_tolerance>`` element specifies two parameters. The first is
the absolute inner tolerance for Gauss-Seidel iterations when performing CMFD
and the second is the relative inner tolerance for Gauss-Seidel iterations
for CMFD calculations.
*Default*: 1.e-10 1.e-5
--------------------
``<ktol>`` Element
--------------------
The ``<ktol>`` element specifies the tolerance on the eigenvalue when performing
CMFD power iteration.
*Default*: 1.e-8
------------------
``<mesh>`` Element
------------------
The CMFD mesh is a structured Cartesian mesh. This element has the following
attributes/sub-elements:
:lower_left:
The lower-left corner of the structured mesh. If only two coordinates are
given, it is assumed that the mesh is an x-y mesh.
:upper_right:
The upper-right corner of the structrued mesh. If only two coordinates are
given, it is assumed that the mesh is an x-y mesh.
:dimension:
The number of mesh cells in each direction.
:width:
The width of mesh cells in each direction.
:energy:
Energy bins [in eV], listed in ascending order (e.g. 0.0 0.625 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:
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.
*Default*: 1.0 1.0 1.0 1.0 1.0 1.0
:map:
An optional acceleration map can be specified to overlay on the coarse
mesh spatial grid. If this option is used, a ``1`` is used for a
non-accelerated region and a ``2`` is used for an accelerated region.
For a simple 4x4 coarse mesh with a 2x2 fuel lattice surrounded by
reflector, the map is:
``1 1 1 1``
``1 2 2 1``
``1 2 2 1``
``1 1 1 1``
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 ``2`` must be used to identify any fission source region.
.. note:: Only two of the following three sub-elements are needed:
``lower_left``, ``upper_right`` and ``width``. Any combination
of two of these will yield the third.
------------------
``<norm>`` Element
------------------
The ``<norm>`` element is used to normalize the CMFD fission source distribution
to a particular value. For example, if a fission source is calculated for a
17 x 17 lattice of pins, the fission source may be normalized to the number of
fission source regions, in this case 289. This is useful when visualizing this
distribution as the average peaking factor will be unity. This parameter will
not impact the calculation.
*Default*: 1.0
---------------------------
``<power_monitor>`` Element
---------------------------
The ``<power_monitor>`` element is used to view the convergence of power
iteration. This option can be turned on with "true" and turned off with "false".
*Default*: false
-------------------------
``<run_adjoint>`` Element
-------------------------
The ``<run_adjoint>`` element can be turned on with "true" to have an adjoint
calculation be performed on the last batch when CMFD is active.
*Default*: false
--------------------
``<shift>`` Element
--------------------
The ``<shift>`` element specifies an optional Wielandt shift parameter for
accelerating power iterations. It is by default very large so the impact of the
shift is effectively zero.
*Default*: 1e6
----------------------
``<spectral>`` Element
----------------------
The ``<spectral>`` element specifies an optional spectral radius that can be set to
accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration
solve.
*Default*: 0.0
------------------
``<stol>`` Element
------------------
The ``<stol>`` element specifies the tolerance on the fission source when performing
CMFD power iteration.
*Default*: 1.e-8
-------------------------
``<tally_reset>`` Element
-------------------------
The ``<tally_reset>`` element contains a list of batch numbers in which CMFD tallies
should be reset.
*Default*: None
----------------------------
``<write_matrices>`` Element
----------------------------
The ``<write_matrices>`` element is used to write the sparse matrices created
when solving CMFD equations. This option can be turned on with "true" and off
with "false".
*Default*: false

1
docs/source/io_formats/index.rst
View file

@ -19,7 +19,6 @@ Input Files
settings
tallies
plots
cmfd
----------
Data Files

10
docs/source/io_formats/settings.rst
View file

@ -361,16 +361,6 @@ or sub-elements:
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
----------------------
``<run_cmfd>`` Element
----------------------
The ``<run_cmfd>`` element indicates whether or not CMFD acceleration should be
turned on or off. This element has no attributes or sub-elements and can be set
to either "false" or "true".
*Default*: false
----------------------
``<run_mode>`` Element
----------------------

20
docs/source/io_formats/statepoint.rst
View file

@ -17,8 +17,6 @@ The current version of the statepoint file format is 17.0.
- **date_and_time** (*char[]*) -- Date and time the summary was
written.
- **path** (*char[]*) -- Path to directory containing input files.
- **cmfd_on** (*int*) -- Flag indicating whether CMFD is on (1) or
off (0).
- **tallies_present** (*int*) -- Flag indicating whether tallies
are present (1) or not (0).
- **source_present** (*int*) -- Flag indicating whether the source
@ -59,18 +57,6 @@ The current version of the statepoint file format is 17.0.
source particle, respectively. Only present when `run_mode` is
'eigenvalue'.
**/cmfd/**
:Datasets: - **indices** (*int[4]*) -- Indices for cmfd mesh (i,j,k,g)
- **k_cmfd** (*double[]*) -- CMFD eigenvalues
- **cmfd_src** (*double[][][][]*) -- CMFD fission source
- **cmfd_entropy** (*double[]*) -- CMFD estimate of Shannon entropy
- **cmfd_balance** (*double[]*) -- RMS of the residual neutron
balance equation on CMFD mesh
- **cmfd_dominance** (*double[]*) -- CMFD estimate of dominance ratio
- **cmfd_srccmp** (*double[]*) -- RMS comparison of difference
between OpenMC and CMFD fission source
**/tallies/**
:Attributes: - **n_tallies** (*int*) -- Number of user-defined tallies.
@ -166,9 +152,3 @@ All values are given in seconds and are measured on the master process.
source sites between processes for load balancing.
- **accumulating tallies** (*double*) -- Time spent communicating
tally results and evaluating their statistics.
- **CMFD** (*double*) -- Time spent evaluating CMFD.
- **CMFD building matrices** (*double*) -- Time spent buliding CMFD
matrices.
- **CMFD solving matrices** (*double*) -- Time spent solving CMFD
matrices.
- **total** (*double*) -- Total time spent in the program.

42
docs/source/pythonapi/base.rst
View file

@ -130,17 +130,6 @@ Constructing Tallies
openmc.Tally
openmc.Tallies
Coarse Mesh Finite Difference Acceleration
------------------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.CMFDMesh
openmc.CMFD
Geometry Plotting
-----------------
@ -183,7 +172,7 @@ The following classes and functions are used for functional expansion reconstruc
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
:template: myclass.rst
openmc.ZernikeRadial
@ -208,3 +197,32 @@ Various classes may be created when performing tally slicing and/or arithmetic:
openmc.arithmetic.AggregateScore
openmc.arithmetic.AggregateNuclide
openmc.arithmetic.AggregateFilter
Coarse Mesh Finite Difference Acceleration
------------------------------------------
CMFD is implemented in OpenMC and allows users to accelerate fission source
convergence during inactive neutron batches. To run CMFD, the CMFDRun class
should be used, and :mod:`from openmc import cmfd` should be included at the
top of the Python input file. Additionally, this class has a dependence on the
C API.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.CMFDMesh
openmc.CMFDRun
At the minimum, a CMFD mesh needs to be specified in order to run CMFD. Once
these properties are set, an OpenMC simulation can be run with CMFD turned on
with the function:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.CMFDRun.run

4
docs/source/usersguide/basics.rst
View file

@ -43,10 +43,6 @@ described below.
This file gives specifications for producing slice or voxel plots of the
geometry.
:ref:`io_cmfd`
This file specifies execution parameters for coarse mesh finite difference
(CMFD) acceleration.
eXtensible Markup Language (XML)
--------------------------------

25
include/openmc/capi.h
View file

@ -124,6 +124,31 @@ extern "C" {
int openmc_zernike_filter_set_params(int32_t index, const double* x,
const double* y, const double* r);
//! Sets the fixed variables that are used for CMFD linear solver
//! \param[in] CSR format index pointer array of loss matrix
//! \param[in] length of indptr
//! \param[in] CSR format index array of loss matrix
//! \param[in] number of non-zero elements in CMFD loss matrix
//! \param[in] dimension n of nxn CMFD loss matrix
//! \param[in] spectral radius of CMFD matrices and tolerances
//! \param[in] indices storing spatial and energy dimensions of CMFD problem
//! \param[in] coremap for problem, storing accelerated regions
extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr,
const int* indices, int n_elements,
int dim, double spectral,
const int* cmfd_indices,
const int* map);
//! Runs a Gauss Seidel linear solver to solve CMFD matrix equations
//! linear solver
//! \param[in] CSR format data array of coefficient matrix
//! \param[in] right hand side vector
//! \param[out] unknown vector
//! \param[in] tolerance on final error
//! \return number of inner iterations required to reach convergence
extern "C" int openmc_run_linsolver(const double* A_data, const double* b,
double* x, double tol);
// Error codes
extern int OPENMC_E_UNASSIGNED;
extern int OPENMC_E_ALLOCATE;

6
include/openmc/constants.h
View file

@ -434,6 +434,12 @@ constexpr int RUN_MODE_PLOTTING {3};
constexpr int RUN_MODE_PARTICLE {4};
constexpr int RUN_MODE_VOLUME {5};
// ============================================================================
// CMFD CONSTANTS
// For non-accelerated regions on coarse mesh overlay
constexpr int CMFD_NOACCEL {-1};
} // namespace openmc
#endif // OPENMC_CONSTANTS_H

2
openmc/__init__.py
View file

@ -1,5 +1,6 @@
from openmc.arithmetic import *
from openmc.cell import *
from openmc.checkvalue import *
from openmc.mesh import *
from openmc.element import *
from openmc.geometry import *
@ -20,7 +21,6 @@ from openmc.trigger import *
from openmc.tally_derivative import *
from openmc.tallies import *
from openmc.mgxs_library import *
from openmc.cmfd import *
from openmc.executor import *
from openmc.statepoint import *
from openmc.summary import *

24
openmc/capi/core.py
View file

@ -20,6 +20,13 @@ class _Bank(Structure):
('delayed_group', c_int)]
# Define input type for numpy arrays that will be passed into C++ functions
# Must be an int or double array, with single dimension that is contiguous
_array_1d_int = np.ctypeslib.ndpointer(dtype=np.int32, ndim=1,
flags='CONTIGUOUS')
_array_1d_dble = np.ctypeslib.ndpointer(dtype=np.double, ndim=1,
flags='CONTIGUOUS')
_dll.openmc_calculate_volumes.restype = c_int
_dll.openmc_calculate_volumes.errcheck = _error_handler
_dll.openmc_finalize.restype = c_int
@ -36,6 +43,10 @@ _dll.openmc_init.errcheck = _error_handler
_dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
_dll.openmc_get_keff.restype = c_int
_dll.openmc_get_keff.errcheck = _error_handler
_init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int,
c_double, _array_1d_int, _array_1d_int]
_dll.openmc_initialize_linsolver.argtypes = _init_linsolver_argtypes
_dll.openmc_initialize_linsolver.restype = None
_dll.openmc_next_batch.argtypes = [POINTER(c_int)]
_dll.openmc_next_batch.restype = c_int
_dll.openmc_next_batch.errcheck = _error_handler
@ -45,6 +56,10 @@ _dll.openmc_run.restype = c_int
_dll.openmc_run.errcheck = _error_handler
_dll.openmc_reset.restype = c_int
_dll.openmc_reset.errcheck = _error_handler
_run_linsolver_argtypes = [_array_1d_dble, _array_1d_dble, _array_1d_dble,
c_double]
_dll.openmc_run_linsolver.argtypes = _run_linsolver_argtypes
_dll.openmc_run_linsolver.restype = c_int
_dll.openmc_source_bank.argtypes = [POINTER(POINTER(_Bank)), POINTER(c_int64)]
_dll.openmc_source_bank.restype = c_int
_dll.openmc_source_bank.errcheck = _error_handler
@ -126,6 +141,7 @@ def find_material(xyz):
else:
return mats[instance.value]
def hard_reset():
"""Reset tallies, timers, and pseudo-random number generator state."""
_dll.openmc_hard_reset()
@ -315,7 +331,7 @@ def statepoint_write(filename=None, write_source=True):
@contextmanager
def run_in_memory(intracomm=None):
def run_in_memory(**kwargs):
"""Provides context manager for calling OpenMC shared library functions.
This function is intended to be used in a 'with' statement and ensures that
@ -331,11 +347,11 @@ def run_in_memory(intracomm=None):
Parameters
----------
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
**kwargs
All keyword arguments are passed to :func:`init`.
"""
init(intracomm=intracomm)
init(**kwargs)
try:
yield
finally:

2675
openmc/cmfd.py
View file

File diff suppressed because it is too large Load diff

18
openmc/model/model.py
View file

@ -25,8 +25,6 @@ class Model(object):
Settings information
tallies : openmc.Tallies, optional
Tallies information
cmfd : openmc.CMFD, optional
CMFD information
plots : openmc.Plots, optional
Plot information
@ -40,19 +38,16 @@ class Model(object):
Settings information
tallies : openmc.Tallies
Tallies information
cmfd : openmc.CMFD
CMFD information
plots : openmc.Plots
Plot information
"""
def __init__(self, geometry=None, materials=None, settings=None,
tallies=None, cmfd=None, plots=None):
tallies=None, plots=None):
self.geometry = openmc.Geometry()
self.materials = openmc.Materials()
self.settings = openmc.Settings()
self.cmfd = cmfd
self.tallies = openmc.Tallies()
self.plots = openmc.Plots()
@ -83,10 +78,6 @@ class Model(object):
def tallies(self):
return self._tallies
@property
def cmfd(self):
return self._cmfd
@property
def plots(self):
return self._plots
@ -121,11 +112,6 @@ class Model(object):
for tally in tallies:
self._tallies.append(tally)
@cmfd.setter
def cmfd(self, cmfd):
check_type('cmfd', cmfd, (openmc.CMFD, type(None)))
self._cmfd = cmfd
@plots.setter
def plots(self, plots):
check_type('plots', plots, Iterable, openmc.Plot)
@ -196,8 +182,6 @@ class Model(object):
if self.tallies:
self.tallies.export_to_xml()
if self.cmfd is not None:
self.cmfd.export_to_xml()
if self.plots:
self.plots.export_to_xml()

18
openmc/settings.py
View file

@ -90,8 +90,6 @@ class Settings(object):
indicates what nuclides the method should be applied to. In its absence,
the method will be applied to all nuclides with 0 K elastic scattering
data present.
run_cmfd : bool
Indicate if coarse mesh finite difference acceleration is to be used
run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'}
The type of calculation to perform (default is 'eigenvalue')
seed : int
@ -182,7 +180,6 @@ class Settings(object):
self._electron_treatment = None
self._photon_transport = None
self._ptables = None
self._run_cmfd = None
self._seed = None
self._survival_biasing = None
@ -279,10 +276,6 @@ class Settings(object):
def photon_transport(self):
return self._photon_transport
@property
def run_cmfd(self):
return self._run_cmfd
@property
def seed(self):
return self._seed
@ -528,11 +521,6 @@ class Settings(object):
cv.check_type('probability tables', ptables, bool)
self._ptables = ptables
@run_cmfd.setter
def run_cmfd(self, run_cmfd):
cv.check_type('run_cmfd', run_cmfd, bool)
self._run_cmfd = run_cmfd
@seed.setter
def seed(self, seed):
cv.check_type('random number generator seed', seed, Integral)
@ -830,11 +818,6 @@ class Settings(object):
element = ET.SubElement(root, "ptables")
element.text = str(self._ptables).lower()
def _create_run_cmfd_subelement(self, root):
if self._run_cmfd is not None:
element = ET.SubElement(root, "run_cmfd")
element.text = str(self._run_cmfd).lower()
def _create_seed_subelement(self, root):
if self._seed is not None:
element = ET.SubElement(root, "seed")
@ -991,7 +974,6 @@ class Settings(object):
self._create_max_order_subelement(root_element)
self._create_photon_transport_subelement(root_element)
self._create_ptables_subelement(root_element)
self._create_run_cmfd_subelement(root_element)
self._create_seed_subelement(root_element)
self._create_survival_biasing_subelement(root_element)
self._create_cutoff_subelement(root_element)

1
schemas.xml
View file

@ -5,6 +5,5 @@
<documentElement localName="settings" uri="src/relaxng/settings.rnc"/>
<documentElement localName="tallies" uri="src/relaxng/tallies.rnc"/>
<documentElement localName="plots" uri="src/relaxng/plots.rnc"/>
<documentElement localName="cmfd" uri="src/relaxng/cmfd.rnc"/>
<documentElement localName="cross_sections" uri="src/relaxng/cross_sections.rnc"/>
</locatingRules>

8
src/api.F90
View file

@ -59,7 +59,6 @@ contains
subroutine free_memory() bind(C)
use bank_header
use cmfd_header
use geometry_header
use material_header
use photon_header
@ -85,6 +84,9 @@ contains
subroutine free_memory_bank() bind(C)
end subroutine free_memory_bank
subroutine free_memory_cmfd() bind(C)
end subroutine free_memory_cmfd
end interface
call free_memory_geometry()
@ -102,13 +104,11 @@ contains
call free_memory_tally_filter()
call free_memory_tally_derivative()
call free_memory_bank()
call free_memory_cmfd()
#ifdef DAGMC
call free_memory_dagmc()
#endif
! Deallocate CMFD
call deallocate_cmfd(cmfd)
end subroutine free_memory
end module openmc_api

943
src/cmfd_data.F90
View file

@ -1,943 +0,0 @@
module cmfd_data
!==============================================================================
! CMFD_DATA -- This module processes the cmfd tally object to generate
! parameters for CMFD calculation.
!==============================================================================
use cmfd_header, only: allocate_cmfd, cmfd, cmfd_coremap, &
cmfd_downscatter, cmfd_tallies, dhat_reset
use constants
use tally_filter_mesh, only: MeshFilter
implicit none
private
public :: set_up_cmfd, neutron_balance
contains
!==============================================================================
! SET_UP_CMFD configures cmfd object for a CMFD eigenvalue calculation
!==============================================================================
subroutine set_up_cmfd()
use constants, only: CMFD_NOACCEL
! Check for core map and set it up
if ((cmfd_coremap) .and. (cmfd%mat_dim == CMFD_NOACCEL)) call set_coremap()
! Calculate all cross sections based on reaction rates from last batch
call compute_xs()
! Compute effective downscatter cross section
if (cmfd_downscatter) call compute_effective_downscatter()
! Check neutron balance
call neutron_balance()
! Calculate dtilde
call compute_dtilde()
! Calculate dhat
call compute_dhat()
end subroutine set_up_cmfd
!===============================================================================
! COMPUTE_XS takes tallies and computes macroscopic cross sections
!===============================================================================
subroutine compute_xs()
use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, &
FILTER_SURFACE, OUT_LEFT, OUT_RIGHT, OUT_BACK, &
OUT_FRONT, OUT_BOTTOM, OUT_TOP, IN_LEFT, IN_RIGHT, &
IN_BACK, IN_FRONT, IN_BOTTOM, IN_TOP, CMFD_NOACCEL, &
ZERO, ONE, TINY_BIT
use error, only: fatal_error
use mesh_header, only: RegularMesh, meshes
use string, only: to_str
use tally_filter_header, only: filters, filter_matches
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for g
integer :: h ! iteration counter for outgoing groups
integer :: l ! iteration counter for tally filters
integer :: ital ! tally object index
integer :: ijk(3) ! indices for mesh cell
integer :: score_index ! index to pull from tally object
integer :: i_filter_mesh ! index for mesh filter
integer :: i_filter_ein ! index for incoming energy filter
integer :: i_filter_eout ! index for outgoing energy filter
integer :: i_filter_legendre ! index for Legendre filter
integer :: i_mesh ! flattend index for mesh
logical :: energy_filters! energy filters present
real(8) :: flux ! temp variable for flux
type(RegularMesh) :: m ! pointer for mesh object
! Extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! Set flux object and source distribution to all zeros
cmfd % flux = ZERO
cmfd % openmc_src = ZERO
! Associate tallies and mesh
associate (t => cmfd_tallies(1) % obj)
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
end associate
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m = meshes(filt % mesh())
end select
! Set mesh widths
cmfd % hxyz(1,:,:,:) = m % width(1) ! set x width
cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width
cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width
cmfd % keff_bal = ZERO
! Begin loop around tallies
TAL: do ital = 1, size(cmfd_tallies)
! Associate tallies and mesh
associate (t => cmfd_tallies(ital) % obj)
if (ital < 3) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
else if (ital == 3) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESHSURFACE))
else if (ital == 4) then
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_legendre = t % filter(t % find_filter(FILTER_LEGENDRE))
end if
! Check for energy filters
energy_filters = (t % find_filter(FILTER_ENERGYIN) > 0)
if (energy_filters) then
i_filter_ein = t % filter(t % find_filter(FILTER_ENERGYIN))
i_filter_eout = t % filter(t % find_filter(FILTER_ENERGYOUT))
end if
! Begin loop around space
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
! Check for active mesh cell
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! Loop around energy groups
OUTGROUP: do h = 1,ng
! Start tally 1
TALLY: if (ital == 1) then
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
! Apply energy in filter
if (energy_filters) then
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins_data(1) - 1) * t % stride(l)
end do
! Get flux
flux = t % results(RESULT_SUM,1,score_index)
cmfd % flux(h,i,j,k) = flux
! Detect zero flux, abort if located
if ((flux - ZERO) < TINY_BIT) then
call fatal_error('Detected zero flux without coremap overlay &
&at: (' // to_str(i) // ',' // to_str(j) // ',' // &
&to_str(k) // ') in group ' // to_str(h))
end if
! Get total rr and convert to total xs
cmfd % totalxs(h,i,j,k) = t % results(RESULT_SUM,2,score_index) / flux
else if (ital == 2) then
! Begin loop to get energy out tallies
INGROUP: do g = 1, ng
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
if (energy_filters) then
! Apply energy in filter
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
! Set energy out bin
call filter_matches(i_filter_eout) % bins_set_data(1, ng - g + 1)
end if
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins_data(1) - 1) * t % stride(l)
end do
! Get scattering
cmfd % scattxs(h,g,i,j,k) = t % results(RESULT_SUM,1,score_index) /&
cmfd % flux(h,i,j,k)
! Get nu-fission
cmfd % nfissxs(h,g,i,j,k) = t % results(RESULT_SUM,2,score_index) /&
cmfd % flux(h,i,j,k)
! Bank source
cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + &
t % results(RESULT_SUM,2,score_index)
cmfd % keff_bal = cmfd % keff_bal + &
t % results(RESULT_SUM,2,score_index) / t % n_realizations
end do INGROUP
else if (ital == 3) then
! Initialize and filter for energy
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set the bin for this mesh cell
i_mesh = m % get_bin_from_indices([ i, j, k ])
call filter_matches(i_filter_mesh) % bins_set_data(1, 12*(i_mesh - 1) + 1)
! Set the energy bin if needed
if (energy_filters) then
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
score_index = 0
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins_data(1) - 1) * t % stride(l)
end do
! Left surface
cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_LEFT - 1))
cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_LEFT - 1))
! Right surface
cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_RIGHT - 1))
cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_RIGHT - 1))
! Back surface
cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_BACK - 1))
cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_BACK - 1))
! Front surface
cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_FRONT - 1))
cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_FRONT - 1))
! Left surface
cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_BOTTOM - 1))
cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_BOTTOM - 1))
! Right surface
cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(IN_TOP - 1))
cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, &
score_index + 1 + ng*(OUT_TOP - 1))
else if (ital == 4) then
! Reset all bins to 1
do l = 1, size(t % filter)
call filter_matches(t % filter(l)) % bins_clear()
call filter_matches(t % filter(l)) % bins_push_back(1)
end do
! Set ijk as mesh indices
ijk = (/ i, j, k /)
! Get bin number for mesh indices
call filter_matches(i_filter_mesh) &
% bins_set_data(1, m % get_bin_from_indices(ijk))
! Apply energy in filter
if (energy_filters) then
call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
end if
! Apply Legendre filter
call filter_matches(i_filter_legendre) % bins_set_data(1, 2)
! Calculate score index from bins
score_index = 1
do l = 1, size(t % filter)
score_index = score_index + (filter_matches(t % filter(l)) &
% bins_data(1) - 1) * t % stride(l)
end do
! Get p1 scatter rr and convert to p1 scatter xs
cmfd % p1scattxs(h,i,j,k) = &
t % results(RESULT_SUM,1,score_index) / &
cmfd % flux(h,i,j,k)
! Calculate diffusion coefficient
cmfd % diffcof(h,i,j,k) = &
ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - &
cmfd % p1scattxs(h,i,j,k)))
end if TALLY
end do OUTGROUP
end do XLOOP
end do YLOOP
end do ZLOOP
end associate
end do TAL
! Normalize openmc source distribution
cmfd % openmc_src = cmfd % openmc_src/sum(cmfd % openmc_src)*cmfd%norm
end subroutine compute_xs
!===============================================================================
! SET_COREMAP is a routine that sets the core mapping information
!===============================================================================
subroutine set_coremap()
use constants, only: CMFD_NOACCEL
integer :: counter=1 ! counter for unique fuel assemblies
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
! Extract spatial indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
! Count how many fuel assemblies exist
cmfd % mat_dim = sum(cmfd % coremap - 1)
! Allocate indexmap
if (.not. allocated(cmfd % indexmap)) &
allocate(cmfd % indexmap(cmfd % mat_dim,3))
! Begin loops over spatial indices
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
XLOOP: do i = 1, nx
! Check for reflector
if (cmfd % coremap(i,j,k) == 1) then
! reset value to CMFD no acceleration constant
cmfd % coremap(i,j,k) = CMFD_NOACCEL
else
! Must be a fuel --> give unique id number
cmfd % coremap(i,j,k) = counter
cmfd % indexmap(counter,1) = i
cmfd % indexmap(counter,2) = j
cmfd % indexmap(counter,3) = k
counter = counter + 1
end if
end do XLOOP
end do YLOOP
end do ZLOOP
end subroutine set_coremap
!===============================================================================
! NEUTRON_BALANCE computes the RMS neutron balance over the CMFD mesh
!===============================================================================
subroutine neutron_balance()
use constants, only: ONE, ZERO, CMFD_NOACCEL, CMFD_NORES
use simulation_header, only: keff, current_batch
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for g
integer :: h ! iteration counter for outgoing groups
integer :: l ! iteration counter for leakage
integer :: cnt ! number of locations to count neutron balance
real(8) :: leakage ! leakage term in neutron balance
real(8) :: interactions ! total number of interactions in balance
real(8) :: scattering ! scattering term in neutron balance
real(8) :: fission ! fission term in neutron balance
real(8) :: res ! residual of neutron balance
real(8) :: rms ! RMS of the residual
! Extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! Allocate res dataspace
if (.not. allocated(cmfd%resnb)) allocate(cmfd%resnb(ng,nx,ny,nz))
! Reset rms and cnt
rms = ZERO
cnt = 0
! Begin loop around space and energy groups
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
XLOOP: do i = 1, nx
GROUPG: do g = 1, ng
! Check for active mesh
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) then
cmfd%resnb(g,i,j,k) = CMFD_NORES
cycle
end if
end if
! Get leakage
leakage = ZERO
LEAK: do l = 1, 3
leakage = leakage + ((cmfd % current(4*l,g,i,j,k) - &
cmfd % current(4*l-1,g,i,j,k))) - &
((cmfd % current(4*l-2,g,i,j,k) - &
cmfd % current(4*l-3,g,i,j,k)))
end do LEAK
! Interactions
interactions = cmfd % totalxs(g,i,j,k) * cmfd % flux(g,i,j,k)
! Get scattering and fission
scattering = ZERO
fission = ZERO
GROUPH: do h = 1, ng
scattering = scattering + cmfd % scattxs(h,g,i,j,k) * &
cmfd % flux(h,i,j,k)
fission = fission + cmfd % nfissxs(h,g,i,j,k) * &
cmfd % flux(h,i,j,k)
end do GROUPH
! Compute residual
res = leakage + interactions - scattering - (ONE/keff)*fission
! Normalize by flux
res = res/cmfd%flux(g,i,j,k)
! Bank res in cmfd object
cmfd%resnb(g,i,j,k) = res
! Take square for RMS calculation
rms = rms + res**2
cnt = cnt + 1
end do GROUPG
end do XLOOP
end do YLOOP
end do ZLOOP
! Calculate RMS and record in vector for this batch
cmfd % balance(current_batch) = sqrt(ONE/dble(cnt)*rms)
end subroutine neutron_balance
!===============================================================================
! COMPUTE_DTILDE precomputes the diffusion coupling coefficient
!===============================================================================
subroutine compute_dtilde()
use constants, only: CMFD_NOACCEL, ZERO_FLUX, TINY_BIT
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: nxyz(3,2) ! single vector containing boundary locations
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: bound(6) ! vector containing indices for boudary check
real(8) :: albedo(6) ! albedo vector with global boundaries
real(8) :: cell_dc ! diffusion coef of current cell
real(8) :: cell_hxyz(3) ! cell dimensions of current ijk cell
real(8) :: neig_dc ! diffusion coefficient of neighbor cell
real(8) :: neig_hxyz(3) ! cell dimensions of neighbor cell
real(8) :: dtilde ! finite difference coupling parameter
real(8) :: ref_albedo ! albedo to reflector
! Get maximum of spatial and group indices
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! Get boundary condition information
albedo = cmfd%albedo
! Loop over group and spatial indices
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
XLOOP: do i = 1, nx
GROUP: do g = 1, ng
! Check for active mesh cell
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) cycle
end if
! Get cell data
cell_dc = cmfd%diffcof(g,i,j,k)
cell_hxyz = cmfd%hxyz(:,i,j,k)
! Setup of vector to identify boundary conditions
bound = (/i,i,j,j,k,k/)
! Begin loop around sides of cell for leakage
LEAK: do l = 1, 6
! Define xyz and +/- indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
shift_idx = -2*mod(l,2) + 1 ! shift neig by -1 or +1
! Check if at a boundary
if (bound(l) == nxyz(xyz_idx,dir_idx)) then
! Compute dtilde with albedo boundary condition
dtilde = (2*cell_dc*(1-albedo(l)))/(4*cell_dc*(1+albedo(l)) + &
(1-albedo(l))*cell_hxyz(xyz_idx))
! Check for zero flux
if (abs(albedo(l) - ZERO_FLUX) < TINY_BIT) dtilde = 2*cell_dc / &
cell_hxyz(xyz_idx)
else ! not a boundary
! Compute neighboring cell indices
neig_idx = (/i,j,k/) ! begin with i,j,k
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! Get neigbor cell data
neig_dc = cmfd%diffcof(g,neig_idx(1),neig_idx(2),neig_idx(3))
neig_hxyz = cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3))
! Check for fuel-reflector interface
if (cmfd_coremap) then
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == &
CMFD_NOACCEL .and. cmfd % coremap(i,j,k) /= CMFD_NOACCEL) then
! Get albedo
ref_albedo = get_reflector_albedo(l,g,i,j,k)
! Compute dtilde
dtilde = (2*cell_dc*(1-ref_albedo))/(4*cell_dc*(1+ &
ref_albedo)+(1-ref_albedo)*cell_hxyz(xyz_idx))
else ! Not next to a reflector or no core map
! Compute dtilde to neighbor cell
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
cell_hxyz(xyz_idx)*neig_dc)
end if
else ! no core map
! Compute dtilde to neighbor cell
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
cell_hxyz(xyz_idx)*neig_dc)
end if
end if
! Record dtilde in cmfd object
cmfd%dtilde(l,g,i,j,k) = dtilde
end do LEAK
end do GROUP
end do XLOOP
end do YLOOP
end do ZLOOP
end subroutine compute_dtilde
!===============================================================================
! COMPUTE_DHAT computes the nonlinear coupling coefficient
!===============================================================================
subroutine compute_dhat()
use constants, only: CMFD_NOACCEL, ZERO
use error, only: write_message
use string, only: to_str
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: nxyz(3,2) ! single vector containing boundary locations
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: bound(6) ! vector containing indices for boudary check
real(8) :: cell_dtilde(6) ! cell dtilde for each face
real(8) :: cell_flux ! flux in current cell
real(8) :: current(12) ! area integrated cell current at each face
real(8) :: net_current ! net current on a face
real(8) :: neig_flux ! flux in neighbor cell
real(8) :: dhat ! dhat equivalence parameter
! Get maximum of spatial and group indices
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! Geting loop over group and spatial indices
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
GROUP: do g = 1,ng
! Check for active mesh cell
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! Get cell data
cell_dtilde = cmfd%dtilde(:,g,i,j,k)
cell_flux = cmfd%flux(g,i,j,k)/product(cmfd%hxyz(:,i,j,k))
current = cmfd%current(:,g,i,j,k)
! Setup of vector to identify boundary conditions
bound = (/i,i,j,j,k,k/)
! Begin loop around sides of cell for leakage
LEAK: do l = 1,6
! Define xyz and +/- indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
! Calculate net current on l face (divided by surf area)
net_current = (current(2*l) - current(2*l-1)) / &
product(cmfd%hxyz(:,i,j,k)) * cmfd%hxyz(xyz_idx,i,j,k)
! Check if at a boundary
if (bound(l) == nxyz(xyz_idx,dir_idx)) then
! Compute dhat
dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) / &
cell_flux
else ! not a boundary
! Compute neighboring cell indices
neig_idx = (/i,j,k/) ! begin with i,j,k
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! Get neigbor flux
neig_flux = cmfd%flux(g,neig_idx(1),neig_idx(2),neig_idx(3)) / &
product(cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3)))
! Check for fuel-reflector interface
if (cmfd_coremap) then
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == &
CMFD_NOACCEL .and. cmfd % coremap(i,j,k) /= CMFD_NOACCEL) then
! compute dhat
dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) /&
cell_flux
else ! not a fuel-reflector interface
! Compute dhat
dhat = (net_current + shift_idx*cell_dtilde(l)* &
(neig_flux - cell_flux))/(neig_flux + cell_flux)
end if
else ! not for fuel-reflector case
! Compute dhat
dhat = (net_current + shift_idx*cell_dtilde(l)* &
(neig_flux - cell_flux))/(neig_flux + cell_flux)
end if
end if
! record dhat in cmfd object
cmfd%dhat(l,g,i,j,k) = dhat
! check for dhat reset
if (dhat_reset) then
cmfd%dhat(l,g,i,j,k) = ZERO
end if
end do LEAK
end do GROUP
end do XLOOP
end do YLOOP
end do ZLOOP
! write that dhats are zero
if (dhat_reset) then
call write_message('Dhats reset to zero.', 8)
end if
end subroutine compute_dhat
!===============================================================================
! GET_REFLECTOR_ALBEDO is a function that calculates the albedo to the reflector
!===============================================================================
function get_reflector_albedo(l, g, i, j, k)
use constants, only: ONE
real(8) :: get_reflector_albedo ! reflector albedo
integer, intent(in) :: i ! iteration counter for x
integer, intent(in) :: j ! iteration counter for y
integer, intent(in) :: k ! iteration counter for z
integer, intent(in) :: g ! iteration counter for groups
integer, intent(in) :: l ! iteration counter for leakages
integer :: shift_idx ! parameter to shift index by +1 or -1
real(8) :: current(12) ! partial currents for all faces of mesh cell
real(8) :: albedo ! the albedo
! Get partial currents from object
current = cmfd%current(:,g,i,j,k)
! Define xyz and +/- indices
shift_idx = -2*mod(l,2) + 1 ! shift neig by -1 or +1
! Calculate albedo
if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. &
(shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then
albedo = ONE
else
albedo = (current(2*l-1)/current(2*l))**(shift_idx)
end if
! Assign to function variable
get_reflector_albedo = albedo
end function get_reflector_albedo
!===============================================================================
! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter
!===============================================================================
subroutine compute_effective_downscatter()
use constants, only: ZERO, CMFD_NOACCEL
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
real(8) :: flux1 ! group 1 volume int flux
real(8) :: flux2 ! group 2 volume int flux
real(8) :: sigt1 ! group 1 total xs
real(8) :: sigt2 ! group 2 total xs
real(8) :: sigs11 ! scattering transfer 1 --> 1
real(8) :: sigs21 ! scattering transfer 2 --> 1
real(8) :: sigs12 ! scattering transfer 1 --> 2
real(8) :: sigs22 ! scattering transfer 2 --> 2
real(8) :: siga1 ! group 1 abs xs
real(8) :: siga2 ! group 2 abs xs
real(8) :: sigs12_eff ! effective downscatter xs
! Extract spatial and energy indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! Return if not two groups
if (ng /= 2) return
! Begin loop around space and energy groups
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
XLOOP: do i = 1, nx
! Check for active mesh
if (allocated(cmfd%coremap)) then
if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) cycle
end if
! Extract cross sections and flux from object
flux1 = cmfd % flux(1,i,j,k)
flux2 = cmfd % flux(2,i,j,k)
sigt1 = cmfd % totalxs(1,i,j,k)
sigt2 = cmfd % totalxs(2,i,j,k)
sigs11 = cmfd % scattxs(1,1,i,j,k)
sigs21 = cmfd % scattxs(2,1,i,j,k)
sigs12 = cmfd % scattxs(1,2,i,j,k)
sigs22 = cmfd % scattxs(2,2,i,j,k)
! Compute absorption xs
siga1 = sigt1 - sigs11 - sigs12
siga2 = sigt2 - sigs22 - sigs21
! Compute effective downscatter xs
sigs12_eff = sigs12 - sigs21*flux2/flux1
! Recompute total cross sections (use effective and no upscattering)
sigt1 = siga1 + sigs11 + sigs12_eff
sigt2 = siga2 + sigs22
! Record total xs
cmfd % totalxs(1,i,j,k) = sigt1
cmfd % totalxs(2,i,j,k) = sigt2
! Record effective downscatter xs
cmfd % scattxs(1,2,i,j,k) = sigs12_eff
! Zero out upscatter cross section
cmfd % scattxs(2,1,i,j,k) = ZERO
end do XLOOP
end do YLOOP
end do ZLOOP
end subroutine compute_effective_downscatter
end module cmfd_data

410
src/cmfd_execute.F90
View file

@ -1,410 +0,0 @@
module cmfd_execute
!==============================================================================
! CMFD_EXECUTE -- This module is the highest level cmfd module that controls the
! cross section generation, diffusion calculation, and source re-weighting
!==============================================================================
use cmfd_header
use settings
use simulation_header
implicit none
private
public :: execute_cmfd, cmfd_init_batch, cmfd_tally_init
#ifdef OPENMC_MPI
interface
subroutine cmfd_broadcast(n, buffer) bind(C)
import C_DOUBLE, C_INT
integer(C_INT), value :: n
real(C_DOUBLE), intent(out) :: buffer
end subroutine
end interface
#endif
contains
!==============================================================================
! EXECUTE_CMFD runs the CMFD calculation
!==============================================================================
subroutine execute_cmfd() bind(C)
use cmfd_data, only: set_up_cmfd
use cmfd_solver, only: cmfd_solver_execute
use error, only: warning, fatal_error
use message_passing, only: master
! CMFD single processor on master
if (master) then
! Start cmfd timer
call time_cmfd % start()
! Create cmfd data from OpenMC tallies
call set_up_cmfd()
! Call solver
call cmfd_solver_execute()
! Save k-effective
cmfd % k_cmfd(current_batch) = cmfd % keff
! check to perform adjoint on last batch
if (current_batch == n_batches .and. cmfd_run_adjoint) then
call cmfd_solver_execute(adjoint=.true.)
end if
end if
! calculate fission source
call calc_fission_source()
! calculate weight factors
call cmfd_reweight(.true.)
! stop cmfd timer
if (master) call time_cmfd % stop()
end subroutine execute_cmfd
!==============================================================================
! CMFD_INIT_BATCH handles cmfd options at the start of every batch
!==============================================================================
subroutine cmfd_init_batch() bind(C)
! Check to activate CMFD diffusion and possible feedback
! this guarantees that when cmfd begins at least one batch of tallies are
! accumulated
if (cmfd_run .and. cmfd_begin == current_batch) then
cmfd_on = .true.
end if
! If this is a restart run and we are just replaying batches leave
if (restart_run .and. current_batch <= restart_batch) return
! Check to reset tallies
if (cmfd_run .and. cmfd_reset % contains(current_batch)) then
call cmfd_tally_reset()
end if
end subroutine cmfd_init_batch
!===============================================================================
! CALC_FISSION_SOURCE calculates the cmfd fission source
!===============================================================================
subroutine calc_fission_source()
use constants, only: CMFD_NOACCEL, ZERO, TWO
use message_passing
use string, only: to_str
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer :: ng ! maximum number of energy groups
integer :: n ! total size
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: idx ! index in vector
real(8) :: hxyz(3) ! cell dimensions of current ijk cell
real(8) :: vol ! volume of cell
real(8),allocatable :: source(:,:,:,:) ! tmp source array for entropy
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
n = ng*nx*ny*nz
! Allocate cmfd source if not already allocated and allocate buffer
if (.not. allocated(cmfd % cmfd_src)) &
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
! Reset cmfd source to 0
cmfd % cmfd_src = ZERO
! Only perform for master
if (master) then
! Loop around indices to map to cmfd object
ZLOOP: do k = 1, nz
YLOOP: do j = 1, ny
XLOOP: do i = 1, nx
GROUP: do g = 1, ng
! Check for core map
if (cmfd_coremap) then
if (cmfd % coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! Get dimensions of cell
hxyz = cmfd % hxyz(:,i,j,k)
! Calculate volume
vol = hxyz(1)*hxyz(2)*hxyz(3)
! Get first index
idx = get_matrix_idx(1,i,j,k,ng,nx,ny)
! Compute fission source
cmfd % cmfd_src(g,i,j,k) = sum(cmfd % nfissxs(:,g,i,j,k) * &
cmfd % phi(idx:idx + (ng - 1)))*vol
end do GROUP
end do XLOOP
end do YLOOP
end do ZLOOP
! Normalize source such that it sums to 1.0
cmfd % cmfd_src = cmfd % cmfd_src/sum(cmfd % cmfd_src)
! Compute entropy
if (entropy_on) then
! Allocate tmp array
if (.not.allocated(source)) allocate(source(ng,nx,ny,nz))
! Initialize the source
source = ZERO
! Compute log
where (cmfd % cmfd_src > ZERO)
source = cmfd % cmfd_src*log(cmfd % cmfd_src)/log(TWO)
end where
! Sum that source
cmfd % entropy(current_batch) = -sum(source)
! Deallocate tmp array
if (allocated(source)) deallocate(source)
end if
! Normalize source so average is 1.0
cmfd % cmfd_src = cmfd % cmfd_src/sum(cmfd % cmfd_src)*cmfd % norm
! Calculate differences between normalized sources
cmfd % src_cmp(current_batch) = sqrt(ONE/cmfd % norm * &
sum((cmfd % cmfd_src - cmfd % openmc_src)**2))
end if
#ifdef OPENMC_MPI
! Broadcast full source to all procs
call cmfd_broadcast(n, cmfd % cmfd_src(1,1,1,1))
#endif
end subroutine calc_fission_source
!===============================================================================
! CMFD_REWEIGHT performs weighting of particles in the source bank
!===============================================================================
subroutine cmfd_reweight(new_weights)
use algorithm, only: binary_search
use bank_header
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use message_passing
use string, only: to_str
logical, intent(in) :: new_weights ! calcualte new weights
integer :: nx ! maximum number of cells in x direction
integer :: ny ! maximum number of cells in y direction
integer :: nz ! maximum number of cells in z direction
integer(C_INT) :: ng ! maximum number of energy groups
integer(8) :: i ! iteration counter
integer :: g ! index for group
integer :: ijk(3) ! spatial bin location
integer :: e_bin ! energy bin of source particle
integer :: mesh_bin ! mesh bin of soruce particle
integer :: n_groups ! number of energy groups
real(8) :: norm ! normalization factor
real(C_DOUBLE) :: xyz(3)
logical(C_BOOL) :: outside ! any source sites outside mesh
logical :: in_mesh ! source site is inside mesh
interface
subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C)
import C_INT, C_DOUBLE, C_BOOL
integer(C_INT), value :: ng
real(C_DOUBLE), intent(in) :: energies
real(C_DOUBLE), intent(out) :: source_counts
logical(C_BOOL), intent(out) :: outside
end subroutine
end interface
! Get maximum of spatial and group indices
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
ng = cmfd % indices(4)
! allocate arrays in cmfd object (can take out later extend to multigroup)
if (.not.allocated(cmfd%sourcecounts)) then
allocate(cmfd%sourcecounts(ng, nx*ny*nz))
cmfd % sourcecounts = 0
end if
if (.not.allocated(cmfd % weightfactors)) then
allocate(cmfd % weightfactors(ng,nx,ny,nz))
cmfd % weightfactors = ONE
end if
! Compute new weight factors
if (new_weights) then
! Set weight factors to a default 1.0
cmfd%weightfactors = ONE
! Count bank sites in mesh and reverse due to egrid structure
call cmfd_populate_sourcecounts(ng + 1, cmfd % egrid(1), &
cmfd % sourcecounts(1,1), outside)
! Check for sites outside of the mesh
if (master .and. outside) then
call fatal_error("Source sites outside of the CMFD mesh!")
end if
! Have master compute weight factors (watch for 0s)
if (master) then
! Calculate normalization factor
norm = sum(cmfd % sourcecounts) / sum(cmfd % cmfd_src)
do mesh_bin = 1, nx*ny*nz
call cmfd_mesh % get_indices_from_bin(mesh_bin, ijk)
do g = 1, ng
if (cmfd % sourcecounts(ng - g + 1, mesh_bin) > ZERO) then
if (cmfd % cmfd_src(g,ijk(1),ijk(2),ijk(3)) > ZERO) then
cmfd % weightfactors(g,ijk(1),ijk(2),ijk(3)) = &
cmfd % cmfd_src(g,ijk(1),ijk(2),ijk(3)) * norm &
/ cmfd % sourcecounts(ng - g + 1, mesh_bin)
end if
end if
end do
end do
end if
if (.not. cmfd_feedback) return
! Broadcast weight factors to all procs
#ifdef OPENMC_MPI
call cmfd_broadcast(ng*nx*ny*nz, cmfd % weightfactors(1,1,1,1))
#endif
end if
! begin loop over source bank
do i = 1, work
! Determine spatial bin
call source_bank_xyz(i, xyz)
call cmfd_mesh % get_indices(xyz, ijk, in_mesh)
! Determine energy bin
n_groups = size(cmfd % egrid) - 1
if (source_bank_E(i) < cmfd % egrid(1)) then
e_bin = 1
if (master) call warning('Source pt below energy grid')
elseif (source_bank_E(i) > cmfd % egrid(n_groups + 1)) then
e_bin = n_groups
if (master) call warning('Source pt above energy grid')
else
e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank_E(i))
end if
! Reverese energy bin (lowest grp is highest energy bin)
e_bin = n_groups - e_bin + 1
! Check for outside of mesh
if (.not. in_mesh) then
call fatal_error('Source site found outside of CMFD mesh')
end if
! Reweight particle
call source_bank_set_wgt(i, source_bank_wgt(i) * &
cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3)))
end do
end subroutine cmfd_reweight
!===============================================================================
! GET_MATRIX_IDX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
function get_matrix_idx(g, i, j, k, ng, nx, ny) result (matidx)
integer :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of cells in x direction
integer, intent(in) :: ny ! maximum number of cells in y direction
integer, intent(in) :: ng ! maximum number of energy groups
! Check if coremap is used
if (cmfd_coremap) then
! Get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! Compute index
matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
end function get_matrix_idx
!===============================================================================
! CMFD_TALLY_INIT
!===============================================================================
subroutine cmfd_tally_init() bind(C)
integer :: i
if (cmfd_run) then
do i = 1, size(cmfd_tallies)
cmfd_tallies(i) % obj % active = .true.
end do
end if
end subroutine cmfd_tally_init
!===============================================================================
! CMFD_TALLY_RESET resets all cmfd tallies
!===============================================================================
subroutine cmfd_tally_reset()
use error, only: write_message
integer :: i ! loop counter
! Print message
call write_message("CMFD tallies reset", 6)
! Reset CMFD tallies
do i = 1, size(cmfd_tallies)
cmfd_tallies(i) % obj % n_realizations = 0
cmfd_tallies(i) % obj % results(:,:,:) = ZERO
end do
end subroutine cmfd_tally_reset
end module cmfd_execute

40
src/cmfd_execute.cpp
View file

@ -1,40 +0,0 @@
#include <algorithm> // for copy
#include <cstdint>
#include <iostream>
#include "xtensor/xarray.hpp"
#include "xtensor/xio.hpp"
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/simulation.h"
#include "openmc/settings.h"
namespace openmc {
extern "C" void
cmfd_populate_sourcecounts(int n_energy, const double* energies,
double* source_counts, bool* outside)
{
// Get source counts in each mesh bin / energy bin
auto& m = model::meshes.at(settings::index_cmfd_mesh);
xt::xarray<double> counts = m->count_sites(simulation::work,
simulation::source_bank.data(), n_energy, energies, outside);
// Copy data from the xarray into the source counts array
std::copy(counts.begin(), counts.end(), source_counts);
}
#ifdef OPENMC_MPI
extern "C" void
cmfd_broadcast(int n, double* buffer)
{
MPI_Bcast(buffer, n, MPI_DOUBLE, 0, mpi::intracomm);
}
#endif
} // namespace openmc

272
src/cmfd_header.F90
View file

@ -1,272 +0,0 @@
module cmfd_header
use, intrinsic :: ISO_C_BINDING
use constants, only: CMFD_NOACCEL, ZERO, ONE
use mesh_header, only: RegularMesh
use set_header, only: SetInt
use tally_header, only: TallyContainer
use timer_header, only: Timer
implicit none
private
public :: allocate_cmfd, deallocate_cmfd
type, public :: cmfd_type
! Indices for problem
integer :: indices(4)
! Albedo boundary condition
real(8) :: albedo(6)
! Core overlay map
integer, allocatable :: coremap(:,:,:)
integer, allocatable :: indexmap(:,:)
integer :: mat_dim = CMFD_NOACCEL
! Energy grid
real(C_DOUBLE), allocatable :: egrid(:)
! Cross sections
real(8), allocatable :: totalxs(:,:,:,:)
real(8), allocatable :: p1scattxs(:,:,:,:)
real(8), allocatable :: scattxs(:,:,:,:,:)
real(8), allocatable :: nfissxs(:,:,:,:,:)
! Diffusion coefficient
real(8), allocatable :: diffcof(:,:,:,:)
! Current
real(8), allocatable :: current(:,:,:,:,:)
! Flux
real(8), allocatable :: flux(:,:,:,:)
! Coupling coefficients and equivalence parameters
real(8), allocatable :: dtilde(:,:,:,:,:)
real(8), allocatable :: dhat(:,:,:,:,:)
! Dimensions of mesh cells ([hu,hv,hw],xloc,yloc,zloc)
real(8), allocatable :: hxyz(:,:,:,:)
! Source distributions
real(C_DOUBLE), allocatable :: cmfd_src(:,:,:,:)
real(C_DOUBLE), allocatable :: openmc_src(:,:,:,:)
! Source sites in each mesh box
real(C_DOUBLE), allocatable :: sourcecounts(:,:)
! Weight adjustment factors
real(8), allocatable :: weightfactors(:,:,:,:)
! Eigenvector/eigenvalue from cmfd run
real(8), allocatable :: phi(:)
real(8) :: keff = ZERO
! Eigenvector/eigenvalue from adjoint run
real(8), allocatable :: adj_phi(:)
real(8) :: adj_keff = ZERO
! Residual for neutron balance
real(8), allocatable :: resnb(:,:,:,:)
! Openmc source normalization factor
real(8) :: norm = ONE
! "Shannon entropy" from cmfd fission source
real(8), allocatable :: entropy(:)
! RMS of neutron balance equations
real(8), allocatable :: balance(:)
! RMS deviation of OpenMC and CMFD normalized source
real(8), allocatable :: src_cmp(:)
! Dominance ratio
real(8), allocatable :: dom(:)
! List of CMFD k
real(8), allocatable :: k_cmfd(:)
! Balance keff
real(8) :: keff_bal
end type cmfd_type
! Main object
type(cmfd_type), public :: cmfd
integer(C_INT), public, bind(C) :: index_cmfd_mesh
type(RegularMesh), public :: cmfd_mesh
! Pointers for different tallies
type(TallyContainer), public, pointer :: cmfd_tallies(:) => null()
! Timing objects
type(Timer), public :: time_cmfd ! timer for whole cmfd calculation
type(Timer), public :: time_cmfdbuild ! timer for matrix build
type(Timer), public :: time_cmfdsolve ! timer for solver
! Flag for active core map
logical, public :: cmfd_coremap = .false.
! Flag to reset dhats to zero
logical, public :: dhat_reset = .false.
! Flag to activate neutronic feedback via source weights
logical, public :: cmfd_feedback = .false.
! Adjoint method type
character(len=10), public :: cmfd_adjoint_type = 'physical'
! Number of incomplete ilu factorization levels
integer, public :: cmfd_ilu_levels = 1
! Batch to begin cmfd
integer, public :: cmfd_begin = 1
! Tally reset list
integer, public :: n_cmfd_resets
type(SetInt), public :: cmfd_reset
! Compute effective downscatter cross section
logical, public :: cmfd_downscatter = .false.
! Convergence monitoring
logical, public :: cmfd_power_monitor = .false.
! Cmfd output
logical, public :: cmfd_write_matrices = .false.
! Run an adjoint calculation (last batch only)
logical, public :: cmfd_run_adjoint = .false.
! CMFD run logicals
logical(C_BOOL), public, bind(C) :: cmfd_on = .false.
! CMFD display info
character(len=25), public :: cmfd_display = 'balance'
! Estimate of spectral radius of CMFD matrices and tolerances
real(8), public :: cmfd_spectral = ZERO
real(8), public :: cmfd_shift = 1.e6
real(8), public :: cmfd_ktol = 1.e-8_8
real(8), public :: cmfd_stol = 1.e-8_8
real(8), public :: cmfd_atoli = 1.e-10_8
real(8), public :: cmfd_rtoli = 1.e-5_8
contains
!==============================================================================
! ALLOCATE_CMFD allocates all data in of cmfd type
!==============================================================================
subroutine allocate_cmfd(this, n_batches)
integer, intent(in) :: n_batches ! number of batches in calc
type(cmfd_type), intent(inout) :: this ! cmfd instance
integer :: nx ! number of mesh cells in x direction
integer :: ny ! number of mesh cells in y direction
integer :: nz ! number of mesh cells in z direction
integer :: ng ! number of energy groups
! Extract spatial and energy indices from object
nx = this % indices(1)
ny = this % indices(2)
nz = this % indices(3)
ng = this % indices(4)
! Allocate flux, cross sections and diffusion coefficient
if (.not. allocated(this % flux)) allocate(this % flux(ng,nx,ny,nz))
if (.not. allocated(this % totalxs)) allocate(this % totalxs(ng,nx,ny,nz))
if (.not. allocated(this % p1scattxs)) allocate(this % p1scattxs(ng,nx,ny,nz))
if (.not. allocated(this % scattxs)) allocate(this % scattxs(ng,ng,nx,ny,nz))
if (.not. allocated(this % nfissxs)) allocate(this % nfissxs(ng,ng,nx,ny,nz))
if (.not. allocated(this % diffcof)) allocate(this % diffcof(ng,nx,ny,nz))
! Allocate dtilde and dhat
if (.not. allocated(this % dtilde)) allocate(this % dtilde(6,ng,nx,ny,nz))
if (.not. allocated(this % dhat)) allocate(this % dhat(6,ng,nx,ny,nz))
! Allocate dimensions for each box (here for general case)
if (.not. allocated(this % hxyz)) allocate(this % hxyz(3,nx,ny,nz))
! Allocate surface currents
if (.not. allocated(this % current)) allocate(this % current(12,ng,nx,ny,nz))
! Allocate source distributions
if (.not. allocated(this % cmfd_src)) allocate(this % cmfd_src(ng,nx,ny,nz))
if (.not. allocated(this % openmc_src)) allocate(this % openmc_src(ng,nx,ny,nz))
! Allocate source weight modification vars
if (.not. allocated(this % sourcecounts)) allocate(this % sourcecounts(ng,nx*ny*nz))
if (.not. allocated(this % weightfactors)) allocate(this % weightfactors(ng,nx,ny,nz))
! Allocate batchwise parameters
if (.not. allocated(this % entropy)) allocate(this % entropy(n_batches))
if (.not. allocated(this % balance)) allocate(this % balance(n_batches))
if (.not. allocated(this % src_cmp)) allocate(this % src_cmp(n_batches))
if (.not. allocated(this % dom)) allocate(this % dom(n_batches))
if (.not. allocated(this % k_cmfd)) allocate(this % k_cmfd(n_batches))
! Set everthing to 0 except weight multiply factors if feedback isnt on
this % flux = ZERO
this % totalxs = ZERO
this % p1scattxs = ZERO
this % scattxs = ZERO
this % nfissxs = ZERO
this % diffcof = ZERO
this % dtilde = ZERO
this % dhat = ZERO
this % hxyz = ZERO
this % current = ZERO
this % cmfd_src = ZERO
this % openmc_src = ZERO
this % sourcecounts = ZERO
this % weightfactors = ONE
this % balance = ZERO
this % src_cmp = ZERO
this % dom = ZERO
this % k_cmfd = ZERO
this % entropy = ZERO
end subroutine allocate_cmfd
!===============================================================================
! DEALLOCATE_CMFD frees all memory of cmfd type
!===============================================================================
subroutine deallocate_cmfd(this)
type(cmfd_type), intent(inout) :: this ! cmfd instance
if (allocated(this % egrid)) deallocate(this % egrid)
if (allocated(this % totalxs)) deallocate(this % totalxs)
if (allocated(this % p1scattxs)) deallocate(this % p1scattxs)
if (allocated(this % scattxs)) deallocate(this % scattxs)
if (allocated(this % nfissxs)) deallocate(this % nfissxs)
if (allocated(this % diffcof)) deallocate(this % diffcof)
if (allocated(this % current)) deallocate(this % current)
if (allocated(this % flux)) deallocate(this % flux)
if (allocated(this % dtilde)) deallocate(this % dtilde)
if (allocated(this % dhat)) deallocate(this % dhat)
if (allocated(this % hxyz)) deallocate(this % hxyz)
if (allocated(this % coremap)) deallocate(this % coremap)
if (allocated(this % indexmap)) deallocate(this % indexmap)
if (allocated(this % phi)) deallocate(this % phi)
if (allocated(this % sourcecounts)) deallocate(this % sourcecounts)
if (allocated(this % weightfactors)) deallocate(this % weightfactors)
if (allocated(this % cmfd_src)) deallocate(this % cmfd_src)
if (allocated(this % openmc_src)) deallocate(this % openmc_src)
if (allocated(this % balance)) deallocate(this % balance)
if (allocated(this % src_cmp)) deallocate(this % src_cmp)
if (allocated(this % dom)) deallocate(this % dom)
if (allocated(this % k_cmfd)) deallocate(this % k_cmfd)
if (allocated(this % entropy)) deallocate(this % entropy)
if (allocated(this % resnb)) deallocate(this % resnb)
end subroutine deallocate_cmfd
end module cmfd_header

495
src/cmfd_input.F90
View file

@ -1,495 +0,0 @@
module cmfd_input
use, intrinsic :: ISO_C_BINDING
use cmfd_header
use mesh_header
use mgxs_interface, only: energy_bins, num_energy_groups
use tally
use tally_header
use timer_header
implicit none
private
public :: configure_cmfd
contains
!===============================================================================
! CONFIGURE_CMFD initializes CMFD parameters
!===============================================================================
subroutine configure_cmfd()
! Read in cmfd input file
call read_cmfd_xml()
! Initialize timers
call time_cmfd % reset()
call time_cmfdbuild % reset()
call time_cmfdsolve % reset()
! Allocate cmfd object
call allocate_cmfd(cmfd, n_batches)
end subroutine configure_cmfd
!===============================================================================
! READ_INPUT reads the CMFD input file and organizes it into a data structure
!===============================================================================
subroutine read_cmfd_xml()
use constants, only: ZERO, ONE
use error, only: fatal_error, warning, write_message
use string, only: to_lower
use xml_interface
use, intrinsic :: ISO_FORTRAN_ENV
integer :: i, g
integer :: ng
integer :: n_params
integer, allocatable :: iarray(:)
integer, allocatable :: int_array(:)
logical :: file_exists ! does cmfd.xml exist?
logical :: found
character(MAX_LINE_LEN) :: filename
real(8) :: gs_tol(2)
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_mesh
! Read cmfd input file
filename = trim(path_input) // "cmfd.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
! CMFD is optional unless it is in on from settings
if (cmfd_run) then
call fatal_error("No CMFD XML file, '" // trim(filename) // "' does not&
& exist!")
end if
return
else
! Tell user
call write_message("Reading CMFD XML file...", 5)
end if
! Parse cmfd.xml file
call doc % load_file(filename)
root = doc % document_element()
! Get pointer to mesh XML node
node_mesh = root % child("mesh")
! Check if mesh is there
if (.not. node_mesh % associated()) then
call fatal_error("No CMFD mesh specified in CMFD XML file.")
end if
! Set spatial dimensions in cmfd object
call get_node_array(node_mesh, "dimension", cmfd % indices(1:3))
! Get number of energy groups
if (check_for_node(node_mesh, "energy")) then
ng = node_word_count(node_mesh, "energy")
if(.not. allocated(cmfd%egrid)) allocate(cmfd%egrid(ng))
call get_node_array(node_mesh, "energy", cmfd%egrid)
cmfd % indices(4) = ng - 1 ! sets energy group dimension
! If using MG mode, check to see if these egrid points at least match
! the MG Data breakpoints
if (.not. run_CE) then
do i = 1, ng
found = .false.
do g = 1, num_energy_groups + 1
if (cmfd % egrid(i) == energy_bins(g)) then
found = .true.
exit
end if
end do
if (.not. found) then
call fatal_error("CMFD energy mesh boundaries must align with&
& boundaries of multi-group data!")
end if
end do
end if
else
if(.not.allocated(cmfd % egrid)) allocate(cmfd % egrid(2))
cmfd % egrid = [ ZERO, energy_max(NEUTRON) ]
cmfd % indices(4) = 1 ! one energy group
end if
! Set global albedo
if (check_for_node(node_mesh, "albedo")) then
call get_node_array(node_mesh, "albedo", cmfd % albedo)
else
cmfd % albedo = [ ONE, ONE, ONE, ONE, ONE, ONE ]
end if
! Get acceleration map
if (check_for_node(node_mesh, "map")) then
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
cmfd % indices(3)))
if (node_word_count(node_mesh, "map") /= &
product(cmfd % indices(1:3))) then
call fatal_error('CMFD coremap not to correct dimensions')
end if
allocate(iarray(node_word_count(node_mesh, "map")))
call get_node_array(node_mesh, "map", iarray)
cmfd % coremap = reshape(iarray,(cmfd % indices(1:3)))
cmfd_coremap = .true.
deallocate(iarray)
end if
! Check for normalization constant
if (check_for_node(root, "norm")) then
call get_node_value(root, "norm", cmfd % norm)
end if
! Set feedback logical
if (check_for_node(root, "feedback")) then
call get_node_value(root, "feedback", cmfd_feedback)
end if
! Set downscatter logical
if (check_for_node(root, "downscatter")) then
call get_node_value(root, "downscatter", cmfd_downscatter)
end if
! Reset dhat parameters
if (check_for_node(root, "dhat_reset")) then
call get_node_value(root, "dhat_reset", dhat_reset)
end if
! Set monitoring
if (check_for_node(root, "power_monitor")) then
call get_node_value(root, "power_monitor", cmfd_power_monitor)
end if
! Output logicals
if (check_for_node(root, "write_matrices")) then
call get_node_value(root, "write_matrices", cmfd_write_matrices)
end if
! Run an adjoint calc
if (check_for_node(root, "run_adjoint")) then
call get_node_value(root, "run_adjoint", cmfd_run_adjoint)
end if
! Batch to begin cmfd
if (check_for_node(root, "begin")) &
call get_node_value(root, "begin", cmfd_begin)
! Check for cmfd tally resets
if (check_for_node(root, "tally_reset")) then
n_cmfd_resets = node_word_count(root, "tally_reset")
else
n_cmfd_resets = 0
end if
if (n_cmfd_resets > 0) then
allocate(int_array(n_cmfd_resets))
call get_node_array(root, "tally_reset", int_array)
do i = 1, n_cmfd_resets
call cmfd_reset % add(int_array(i))
end do
deallocate(int_array)
end if
! Get display
if (check_for_node(root, "display")) &
call get_node_value(root, "display", cmfd_display)
! Read in spectral radius estimate and tolerances
if (check_for_node(root, "spectral")) &
call get_node_value(root, "spectral", cmfd_spectral)
if (check_for_node(root, "shift")) &
call get_node_value(root, "shift", cmfd_shift)
if (check_for_node(root, "ktol")) &
call get_node_value(root, "ktol", cmfd_ktol)
if (check_for_node(root, "stol")) &
call get_node_value(root, "stol", cmfd_stol)
if (check_for_node(root, "gauss_seidel_tolerance")) then
n_params = node_word_count(root, "gauss_seidel_tolerance")
if (n_params /= 2) then
call fatal_error('Gauss Seidel tolerance is not 2 parameters &
&(absolute, relative).')
end if
call get_node_array(root, "gauss_seidel_tolerance", gs_tol)
cmfd_atoli = gs_tol(1)
cmfd_rtoli = gs_tol(2)
end if
! Create tally objects
call create_cmfd_tally(root)
! Close CMFD XML file
call doc % clear()
end subroutine read_cmfd_xml
!===============================================================================
! CREATE_CMFD_TALLY creates the tally object for OpenMC to process for CMFD
! accleration.
! There are 3 tally types:
! 1: Only an energy in filter-> flux,total,p1 scatter
! 2: Energy in and energy out filter-> nu-scatter,nu-fission
! 3: Mesh current
!===============================================================================
subroutine create_cmfd_tally(root)
use constants, only: MAX_LINE_LEN
use error, only: fatal_error, warning
use mesh_header
use string
use tally, only: openmc_tally_allocate
use tally_header, only: openmc_extend_tallies
use tally_filter_header
use tally_filter
use xml_interface
type(XMLNode), intent(in) :: root ! XML root element
logical :: energy_filters
integer :: i ! loop counter
integer :: n ! size of arrays in mesh specification
integer(C_INT32_T) :: ng ! number of energy groups (default 1)
integer :: n_filter ! number of filters
integer :: i_start, i_end
integer :: i_filt_start, i_filt_end
integer(C_INT32_T), allocatable :: filter_indices(:)
integer(C_INT) :: err
integer :: i_filt ! index in filters array
integer :: filt_id
integer :: tally_id
real(C_DOUBLE), allocatable :: energies(:)
type(XMLNode) :: node_mesh
! Read CMFD mesh
call read_meshes(root % ptr)
! Get index of cmfd mesh and set ID
i_start = n_meshes() - 1
err = openmc_mesh_set_id(i_start, i_start)
! Save reference to CMFD mesh
index_cmfd_mesh = i_start
cmfd_mesh = meshes(i_start)
! Get pointer to mesh XML node
node_mesh = root % child("mesh")
! Determine number of filters
energy_filters = check_for_node(node_mesh, "energy")
n = merge(5, 3, energy_filters)
! Extend filters array so we can add CMFD filters
err = openmc_extend_filters(n, i_filt_start, i_filt_end)
! Set up mesh filter
i_filt = i_filt_start
err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_mesh_filter_set_mesh(i_filt, i_start)
if (energy_filters) then
! Read and set incoming energy mesh filter
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'energy' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
! Get energies and set bins
ng = node_word_count(node_mesh, "energy")
allocate(energies(ng))
call get_node_array(node_mesh, "energy", energies)
err = openmc_energy_filter_set_bins(i_filt, ng, energies)
! Read and set outgoing energy mesh filter
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'energyout' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_energy_filter_set_bins(i_filt, ng, energies)
end if
! Duplicate the mesh filter for the mesh current tally since other
! tallies use this filter and we need to change the dimension
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'meshsurface' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_meshsurface_filter_set_mesh(i_filt, i_start)
! Add in legendre filter for the P1 tally
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'legendre' // C_NULL_CHAR)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_legendre_filter_set_order(i_filt, 1)
! Initialize filters
do i = i_filt_start, i_filt_end
call filters(i) % obj % initialize()
end do
! Allocate tallies
err = openmc_extend_tallies(4, i_start, i_end)
cmfd_tallies => tallies(i_start:i_end)
! Begin loop around tallies
do i = 1, size(cmfd_tallies)
! Allocate tally
err = openmc_tally_allocate(i_start + i - 1, C_CHAR_'generic' // C_NULL_CHAR)
call openmc_get_tally_next_id(tally_id)
err = openmc_tally_set_id(i_start + i - 1, tally_id)
! Point t to tally variable
associate (t => cmfd_tallies(i) % obj)
! Set the incoming energy mesh filter index in the tally find_filter
! array
n_filter = 1
if (energy_filters) then
n_filter = n_filter + 1
end if
! Set number of nucilde bins
allocate(t % nuclide_bins(1))
t % nuclide_bins(1) = -1
t % n_nuclide_bins = 1
! Record tally id which is equivalent to loop number
t % id = i_start + i - 1
if (i == 1) then
! Set name
t % name = "CMFD flux, total"
! Set tally estimator to analog
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Allocate and set filters
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_start
if (energy_filters) then
filter_indices(2) = i_filt_start + 1
end if
err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices)
deallocate(filter_indices)
! Allocate scoring bins
allocate(t % score_bins(2))
t % n_score_bins = 2
! Set macro_bins
t % score_bins(1) = SCORE_FLUX
t % score_bins(2) = SCORE_TOTAL
else if (i == 2) then
! Set name
t % name = "CMFD neutron production"
! Set tally estimator to analog
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Set the incoming energy mesh filter index in the tally find_filter
! array
if (energy_filters) then
n_filter = n_filter + 1
end if
! Allocate and set indices in filters array
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_start
if (energy_filters) then
filter_indices(2) = i_filt_start + 1
filter_indices(3) = i_filt_start + 2
end if
err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices)
deallocate(filter_indices)
! Allocate macro reactions
allocate(t % score_bins(2))
t % n_score_bins = 2
! Set macro_bins
t % score_bins(1) = SCORE_NU_SCATTER
t % score_bins(2) = SCORE_NU_FISSION
else if (i == 3) then
! Set name
t % name = "CMFD surface currents"
! Set tally estimator to analog
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Allocate and set filters
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_end - 1
if (energy_filters) then
filter_indices(2) = i_filt_start + 1
end if
err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices)
deallocate(filter_indices)
! Allocate macro reactions
allocate(t % score_bins(1))
t % n_score_bins = 1
! Set macro bins
t % score_bins(1) = SCORE_CURRENT
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'mesh-surface' // C_NULL_CHAR)
else if (i == 4) then
! Set name
t % name = "CMFD P1 scatter"
! Set tally estimator to analog
err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR)
! Set tally type to volume
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Allocate and set filters
n_filter = 2
if (energy_filters) then
n_filter = n_filter + 1
end if
allocate(filter_indices(n_filter))
filter_indices(1) = i_filt_start
filter_indices(2) = i_filt_end
if (energy_filters) then
filter_indices(3) = i_filt_start + 1
end if
err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices)
deallocate(filter_indices)
! Allocate scoring bins
allocate(t % score_bins(1))
t % n_score_bins = 1
! Set macro_bins
t % score_bins(1) = SCORE_SCATTER
end if
! Make CMFD tallies active from the start
t % active = .true.
end associate
end do
end subroutine create_cmfd_tally
end module cmfd_input

435
src/cmfd_loss_operator.F90
View file

@ -1,435 +0,0 @@
module cmfd_loss_operator
use constants, only: CMFD_NOACCEL, ZERO
use cmfd_header, only: cmfd, cmfd_coremap
use matrix_header, only: Matrix
implicit none
private
public :: init_loss_matrix, build_loss_matrix
contains
!===============================================================================
! INIT_LOSS_MATRIX preallocates loss matrix and initializes it
!===============================================================================
subroutine init_loss_matrix(loss_matrix)
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: n ! total length of matrix
integer :: nnz ! number of nonzeros in matrix
integer :: n_i ! number of interior cells
integer :: n_c ! number of corner cells
integer :: n_s ! number of side cells
integer :: n_e ! number of edge cells
integer :: nz_c ! number of non-zero corner cells
integer :: nz_e ! number of non-zero edge cells
integer :: nz_s ! number of non-zero side cells
integer :: nz_i ! number of non-zero interior cells
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Calculate dimensions of matrix
if (cmfd_coremap) then
n = cmfd % mat_dim * ng
else
n = nx*ny*nz*ng
end if
! Calculate number of nonzeros, if core map -> need to determine manually
if (cmfd_coremap) then
nnz = preallocate_loss_matrix(nx, ny, nz, ng, n)
else ! structured Cartesian grid
n_c = 8 ! define # of corners
n_e = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges
n_s = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) &
+ 2*(ny - 2)*(nz - 2) ! define # of sides
n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors
nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners
nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges
nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides
nz_i = ng*n_i*(7 + ng - 1) ! define # nonzero interiors
nnz = nz_c + nz_e + nz_s + nz_i
end if
! Configure loss matrix
call loss_matrix % create(n, nnz)
end subroutine init_loss_matrix
!===============================================================================
! PREALLOCATE_LOSS_MATRIX manually preallocates the loss matrix
!===============================================================================
function preallocate_loss_matrix(nx, ny, nz, ng, n) result(nnz)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
integer, intent(in) :: n ! total length of matrix
integer :: nnz ! number of nonzeros
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
integer :: irow ! row counter
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: nxyz(3,2) ! single vector containing bound. locations
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
! Reset number of nonzeros to 0
nnz = 0
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! Begin loop around local rows
ROWS: do irow = 1, n
! Set a nonzero for diagonal
nnz = nnz + 1
! Get location indices
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! Begin loop over leakages
LEAK: do l = 1,6
! Define (x,y,z) and (-,+) indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
! Calculate spatial indices of neighbor
neig_idx = (/i,j,k/) ! begin with i,j,k
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! Check for coremap
if (cmfd_coremap) then
! Check for neighbor that is non-acceleartred
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= &
CMFD_NOACCEL) then
! Get neighbor matrix index
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
neig_idx(3), neig_mat_idx, ng, nx, ny)
! Record nonzero
nnz = nnz + 1
end if
else
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx, ng, nx, ny)
! Record nonzero
nnz = nnz + 1
end if
end if
end do LEAK
! Begin loop over off diagonal in-scattering
SCATTR: do h = 1, ng
! Cycle though if h=g, it was already banked in removal xs
if (h == g) cycle
! Get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx, ng, nx, ny)
! Record nonzero
nnz = nnz + 1
end do SCATTR
end do ROWS
end function preallocate_loss_matrix
!===============================================================================
! BUILD_LOSS_MATRIX creates the matrix representing loss of neutrons
!===============================================================================
subroutine build_loss_matrix(loss_matrix, adjoint)
type(Matrix), intent(inout) :: loss_matrix ! cmfd loss matrix
logical, intent(in), optional :: adjoint ! set up the adjoint
integer :: nxyz(3,2) ! single vector containing bound. locations
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: l ! iteration counter for leakages
integer :: h ! energy group when doing scattering
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: neig_mat_idx ! matrix index of neighbor cell
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
integer :: bound(6) ! vector for comparing when looking for bound
integer :: xyz_idx ! index for determining if x,y or z leakage
integer :: dir_idx ! index for determining - or + face of cell
integer :: neig_idx(3) ! spatial indices of neighbour
integer :: shift_idx ! parameter to shift index by +1 or -1
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint calculation?
real(8) :: totxs ! total macro cross section
real(8) :: scattxsgg ! scattering macro cross section g-->g
real(8) :: scattxshg ! scattering macro cross section h-->g
real(8) :: dtilde(6) ! finite difference coupling parameter
real(8) :: dhat(6) ! nonlinear coupling parameter
real(8) :: hxyz(3) ! cell lengths in each direction
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
! Check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Create single vector of these indices for boundary calculation
nxyz(1,:) = (/1,nx/)
nxyz(2,:) = (/1,ny/)
nxyz(3,:) = (/1,nz/)
! Begin iteration loops
ROWS: do irow = 1, loss_matrix % n
! Set up a new row in matrix
call loss_matrix % new_row()
! Get indices for that row
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! Retrieve cell data
totxs = cmfd%totalxs(g,i,j,k)
scattxsgg = cmfd%scattxs(g,g,i,j,k)
dtilde = cmfd%dtilde(:,g,i,j,k)
hxyz = cmfd%hxyz(:,i,j,k)
! Check and get dhat
if (allocated(cmfd%dhat)) then
dhat = cmfd%dhat(:,g,i,j,k)
else
dhat = ZERO
end if
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! Begin loop over leakages
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
LEAK: do l = 1,6
! Define (x,y,z) and (-,+) indices
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
dir_idx = 2 - mod(l,2) ! -=1, +=2
! Calculate spatial indices of neighbor
neig_idx = (/i,j,k/) ! begin with i,j,k
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! Check for core map
if (cmfd_coremap) then
! Check that neighbor is not reflector
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= &
CMFD_NOACCEL) then
! Compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx, ng, nx, ny)
! Compute value and record to bank
val = jn/hxyz(xyz_idx)
! Record value in matrix
call loss_matrix % add_value(neig_mat_idx, val)
end if
else
! Compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! Get neighbor matrix index
call indices_to_matrix(g, neig_idx(1), neig_idx(2), neig_idx(3), &
neig_mat_idx, ng, nx, ny)
! Compute value and record to bank
val = jn/hxyz(xyz_idx)
! Record value in matrix
call loss_matrix % add_value(neig_mat_idx, val)
end if
end if
! Compute leakage coefficient for target
jo(l) = shift_idx*dtilde(l) + dhat(l)
end do LEAK
! Calculate net leakage coefficient for target
jnet = (jo(2) - jo(1))/hxyz(1) + (jo(4) - jo(3))/hxyz(2) + &
(jo(6) - jo(5))/hxyz(3)
! Calculate loss of neutrons
val = jnet + totxs - scattxsgg
! Record diagonal term
call loss_matrix % add_value(irow, val)
! Begin loop over off diagonal in-scattering
SCATTR: do h = 1, ng
! Cycle though if h=g, value already banked in removal xs
if (h == g) cycle
! Get neighbor matrix index
call indices_to_matrix(h, i, j, k, scatt_mat_idx, ng, nx, ny)
! Check for adjoint
if (adjoint_calc) then
! Get scattering macro xs, transposed!
scattxshg = cmfd%scattxs(g, h, i, j, k)
else
! Get scattering macro xs
scattxshg = cmfd%scattxs(h, g, i, j, k)
end if
! Negate the scattering xs
val = -scattxshg
! Record value in matrix
call loss_matrix % add_value(scatt_mat_idx, val)
end do SCATTR
end do ROWS
! CSR requires n+1 row
call loss_matrix % new_row()
end subroutine build_loss_matrix
!===============================================================================
! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
! Check if coremap is used
if (cmfd_coremap) then
! Get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! Compute index
matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
end subroutine indices_to_matrix
!===============================================================================
! MATRIX_TO_INDICES converts a matrix index to spatial and group indices
!===============================================================================
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
! Check for core map
if (cmfd_coremap) then
! Get indices from indexmap
g = mod(irow-1, ng) + 1
i = cmfd % indexmap((irow-1)/ng+1,1)
j = cmfd % indexmap((irow-1)/ng+1,2)
k = cmfd % indexmap((irow-1)/ng+1,3)
else
! Compute indices
g = mod(irow-1, ng) + 1
i = mod(irow-1, ng*nx)/ng + 1
j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1
end if
end subroutine matrix_to_indices
end module cmfd_loss_operator

199
src/cmfd_prod_operator.F90
View file

@ -1,199 +0,0 @@
module cmfd_prod_operator
use constants, only: CMFD_NOACCEL
use cmfd_header, only: cmfd, cmfd_coremap
use matrix_header, only: Matrix
implicit none
private
public :: init_prod_matrix, build_prod_matrix
contains
!==============================================================================
! INIT_PROD_MATRIX preallocates prod matrix and initializes it
!==============================================================================
subroutine init_prod_matrix(prod_matrix)
type(Matrix), intent(inout) :: prod_matrix ! production matrix
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: n ! total length of matrix
integer :: nnz ! number of nonzeros in matrix
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! Calculate dimensions and number of nonzeros in matrix
if (cmfd_coremap) then
n = cmfd % mat_dim * ng
else
n = nx*ny*nz*ng
end if
nnz = n * ng
! Configure prod matrix
call prod_matrix % create(n, nnz)
end subroutine init_prod_matrix
!===============================================================================
! BUILD_PROD_MATRIX creates the matrix representing production of neutrons
!===============================================================================
subroutine build_prod_matrix(prod_matrix, adjoint)
type(Matrix), intent(inout) :: prod_matrix ! production matrix
logical, intent(in), optional :: adjoint ! adjoint calculation logical
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: g ! iteration counter for groups
integer :: h ! energy group when doing scattering
integer :: nx ! maximum number of x cells
integer :: ny ! maximum number of y cells
integer :: nz ! maximum number of z cells
integer :: ng ! maximum number of groups
integer :: hmat_idx ! index in matrix for energy group h
integer :: irow ! iteration counter over row
logical :: adjoint_calc ! is this a physical adjoint?
real(8) :: nfissxs ! nufission cross section h-->g
real(8) :: val ! temporary variable for nfissxs
! get maximum number of cells in each direction
nx = cmfd%indices(1)
ny = cmfd%indices(2)
nz = cmfd%indices(3)
ng = cmfd%indices(4)
! check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! begin iteration loops
ROWS: do irow = 1, prod_matrix % n
! add a new row to matrix
call prod_matrix % new_row()
! get indices for that row
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! check if not including reflector
if (cmfd_coremap) then
! check if at a reflector
if (cmfd % coremap(i,j,k) == CMFD_NOACCEL) then
cycle
end if
end if
! loop around all other groups
NFISS: do h = 1, ng
! get matrix column location
call indices_to_matrix(h, i, j, k, hmat_idx, ng, nx, ny)
! check for adjoint and bank val
if (adjoint_calc) then
! get nu-fission cross section from cell
nfissxs = cmfd%nfissxs(g,h,i,j,k)
else
! get nu-fission cross section from cell
nfissxs = cmfd%nfissxs(h,g,i,j,k)
end if
! set as value to be recorded
val = nfissxs
! record value in matrix
call prod_matrix % add_value(hmat_idx, val)
end do NFISS
end do ROWS
! CSR requires n+1 row
call prod_matrix % new_row()
end subroutine build_prod_matrix
!===============================================================================
! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix
!===============================================================================
subroutine indices_to_matrix(g, i, j, k, matidx, ng, nx, ny)
integer, intent(out) :: matidx ! the index location in matrix
integer, intent(in) :: i ! current x index
integer, intent(in) :: j ! current y index
integer, intent(in) :: k ! current z index
integer, intent(in) :: g ! current group index
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: ng ! maximum number of groups
! check if coremap is used
if (cmfd_coremap) then
! get idx from core map
matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g)
else
! compute index
matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
end subroutine indices_to_matrix
!===============================================================================
! MATRIX_TO_INDICES converts matrix index to spatial and group indices
!===============================================================================
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
! check for core map
if (cmfd_coremap) then
! get indices from indexmap
g = mod(irow-1, ng) + 1
i = cmfd % indexmap((irow-1)/ng+1,1)
j = cmfd % indexmap((irow-1)/ng+1,2)
k = cmfd % indexmap((irow-1)/ng+1,3)
else
! compute indices
g = mod(irow-1, ng) + 1
i = mod(irow-1, ng*nx)/ng + 1
j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1
end if
end subroutine matrix_to_indices
end module cmfd_prod_operator

809
src/cmfd_solver.F90
View file

@ -1,809 +0,0 @@
module cmfd_solver
! This module contains routines to execute the power iteration solver
use constants, only: MAX_LINE_LEN
use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
use matrix_header, only: Matrix
use vector_header, only: Vector
implicit none
private
public :: cmfd_solver_execute
real(8) :: k_n ! New k-eigenvalue
real(8) :: k_o ! Old k-eigenvalue
real(8) :: k_s ! Shift of eigenvalue
real(8) :: k_ln ! New shifted eigenvalue
real(8) :: k_lo ! Old shifted eigenvalue
real(8) :: norm_n ! Current norm of source vector
real(8) :: norm_o ! Old norm of source vector
real(8) :: kerr ! Error in keff
real(8) :: serr ! Error in source
real(8) :: ktol ! Tolerance on keff
real(8) :: stol ! Tolerance on source
logical :: adjoint_calc ! Run an adjoint calculation
type(Matrix) :: loss ! Cmfd loss matrix
type(Matrix) :: prod ! Cmfd prod matrix
type(Vector) :: phi_n ! New flux vector
type(Vector) :: phi_o ! Old flux vector
type(Vector) :: s_n ! New source vector
type(Vector) :: s_o ! Old flux vector
type(Vector) :: serr_v ! Error in source
! CMFD linear solver interface
abstract interface
subroutine linsolve(A, b, x, tol, i)
import :: Matrix
import :: Vector
type(Matrix), intent(inout) :: A
type(Vector), intent(inout) :: b
type(Vector), intent(inout) :: x
real(8), intent(in) :: tol
integer, intent(out) :: i
end subroutine linsolve
end interface
contains
!===============================================================================
! CMFD_SOLVER_EXECUTE sets up and runs power iteration solver for CMFD
!===============================================================================
subroutine cmfd_solver_execute(adjoint)
use cmfd_header, only: cmfd_adjoint_type, time_cmfdbuild, time_cmfdsolve
logical, optional, intent(in) :: adjoint ! adjoint calc
logical :: physical_adjoint = .false.
! Check for adjoint execution
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
! Check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
! Start timer for build
call time_cmfdbuild % start()
! Initialize matrices and vectors
call init_data(physical_adjoint)
! Check for mathematical adjoint calculation
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
! Stop timer for build
call time_cmfdbuild % stop()
! Begin power iteration
call time_cmfdsolve % start()
call execute_power_iter()
call time_cmfdsolve % stop()
! Extract results
call extract_results()
! Deallocate data
call finalize()
end subroutine cmfd_solver_execute
!===============================================================================
! INIT_DATA allocates matrices and vectors for CMFD solution
!===============================================================================
subroutine init_data(adjoint)
use constants, only: ONE, ZERO
use cmfd_header, only: cmfd_shift, cmfd_ktol, cmfd_stol, cmfd_write_matrices
use simulation_header, only: keff
logical, intent(in) :: adjoint
integer :: n ! problem size
real(8) :: guess ! initial guess
real(8) :: dw ! eigenvalue shift
! Set up matrices
call init_loss_matrix(loss)
call init_prod_matrix(prod)
! Get problem size
n = loss % n
! Set up flux vectors
call phi_n % create(n)
call phi_o % create(n)
! Set up source vectors
call s_n % create(n)
call s_o % create(n)
call serr_v % create(n)
! Set initial guess
guess = ONE
phi_n % val = guess
phi_o % val = guess
k_n = keff
k_o = k_n
dw = cmfd_shift
k_s = k_o + dw
k_ln = ONE/(ONE/k_n - ONE/k_s)
k_lo = k_ln
! Fill in loss matrix
call build_loss_matrix(loss, adjoint=adjoint)
! Fill in production matrix
call build_prod_matrix(prod, adjoint=adjoint)
! Finalize setup of CSR matrices
call loss % assemble()
call prod % assemble()
if (cmfd_write_matrices) then
if (adjoint) then
call loss % write('adj_loss.dat')
call prod % write('adj_prod.dat')
else
call loss % write('loss.dat')
call prod % write('prod.dat')
end if
end if
! Set norms to 0
norm_n = ZERO
norm_o = ZERO
! Set tolerances
ktol = cmfd_ktol
stol = cmfd_stol
end subroutine init_data
!===============================================================================
! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem
!===============================================================================
subroutine compute_adjoint()
use error, only: fatal_error
use cmfd_header, only: cmfd_write_matrices
! Transpose matrices
loss = loss % transpose()
prod = prod % transpose()
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call loss % write('adj_loss.dat')
call prod % write('adj_prod.dat')
end if
end subroutine compute_adjoint
!===============================================================================
! EXECUTE_POWER_ITER is the main power iteration routine
! for the cmfd calculation
!===============================================================================
subroutine execute_power_iter()
use constants, only: ONE
use error, only: fatal_error
use cmfd_header, only: cmfd, cmfd_atoli, cmfd_rtoli
integer :: i ! iteration counter
integer :: innerits ! # of inner iterations
integer :: totalits ! total number of inners
logical :: iconv ! did the problem converged
real(8) :: atoli ! absolute minimum tolerance
real(8) :: rtoli ! relative tolerance based on source conv
real(8) :: toli ! the current tolerance of inners
! Reset convergence flag
iconv = .false.
! Set up tolerances
atoli = cmfd_atoli
rtoli = cmfd_rtoli
toli = rtoli*100._8
! Perform shift
call wielandt_shift()
totalits = 0
! Begin power iteration
do i = 1, 10000
! Check if reached iteration 10000
if (i == 10000) then
call fatal_error('Reached maximum iterations in CMFD power iteration &
&solver.')
end if
! Compute source vector
call prod % vector_multiply(phi_o, s_o)
! Normalize source vector
s_o % val = s_o % val / k_lo
! Compute new flux vector
select case(cmfd % indices(4))
case(1)
call cmfd_linsolver_1g(loss, s_o, phi_n, toli, innerits)
case(2)
call cmfd_linsolver_2g(loss, s_o, phi_n, toli, innerits)
case default
call cmfd_linsolver_ng(loss, s_o, phi_n, toli, innerits)
end select
! Compute new source vector
call prod % vector_multiply(phi_n, s_n)
! Compute new shifted eigenvalue
k_ln = sum(s_n % val) / sum(s_o % val)
! Compute new eigenvalue
k_n = ONE/(ONE/k_ln + ONE/k_s)
! Renormalize the old source
s_o % val = s_o % val * k_lo
! Check convergence
call convergence(i, innerits, iconv)
totalits = totalits + innerits
! Break loop if converged
if (iconv) exit
! Record old values
phi_o % val = phi_n % val
k_o = k_n
k_lo = k_ln
norm_o = norm_n
! Get new tolerance for inners
toli = max(atoli, rtoli*serr)
end do
end subroutine execute_power_iter
!===============================================================================
! WIELANDT SHIFT
!===============================================================================
subroutine wielandt_shift()
use constants, only: ONE
integer :: irow ! row counter
integer :: icol ! col counter
integer :: jcol ! current col index in prod matrix
! perform subtraction
jcol = 1
ROWS: do irow = 1, loss % n
COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1
if (loss % get_col(icol) == prod % get_col(jcol) .and. &
jcol < prod % get_row(irow + 1)) then
loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol)
jcol = jcol + 1
end if
end do COLS
end do ROWS
end subroutine wielandt_shift
!===============================================================================
! CONVERGENCE checks the convergence of the CMFD problem
!===============================================================================
subroutine convergence(iter, innerits, iconv)
use, intrinsic :: ISO_FORTRAN_ENV
use constants, only: ONE, ZERO
use cmfd_header, only: cmfd_power_monitor
use message_passing, only: master
integer, intent(in) :: iter ! outer iteration number
integer, intent(in) :: innerits ! inner iteration nubmer
logical, intent(out) :: iconv ! convergence logical
! Reset convergence flag
iconv = .false.
! Calculate error in keff
kerr = abs(k_o - k_n)/k_n
! Calculate max error in source
where (s_n % val > ZERO)
serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2
end where
serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val))
! Check for convergence
if(kerr < ktol .and. serr < stol) iconv = .true.
! Save the L2 norm of the source
norm_n = serr
! Print out to user
if (cmfd_power_monitor .and. master) then
write(OUTPUT_UNIT,FMT='(I0,":",T10,"k-eff: ",F0.8,T30,"k-error: ", &
&1PE12.5,T55, "src-error: ",1PE12.5,T80,I0)') iter, k_n, kerr, &
serr, innerits
end if
end subroutine convergence
!===============================================================================
! CMFD_LINSOLVER_1g solves the CMFD linear system
!===============================================================================
subroutine cmfd_linsolver_1g(A, b, x, tol, its)
use constants, only: ONE, ZERO
use error, only: fatal_error
use cmfd_header, only: cmfd, cmfd_spectral
type(Matrix), intent(inout) :: A ! coefficient matrix
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! unknown vector
real(8), intent(in) :: tol ! tolerance on final error
integer, intent(out) :: its ! number of inner iterations
integer :: g ! group index
integer :: i ! loop counter for x
integer :: j ! loop counter for y
integer :: k ! loop counter for z
integer :: n ! total size of vector
integer :: nx ! maximum dimension in x direction
integer :: ny ! maximum dimension in y direction
integer :: nz ! maximum dimension in z direction
integer :: ng ! number of energy groups
integer :: igs ! Gauss-Seidel iteration counter
integer :: irb ! Red/Black iteration switch
integer :: irow ! row iteration
integer :: icol ! iteration counter over columns
integer :: didx ! index for diagonal component
logical :: found ! did we find col
real(8) :: tmp1 ! temporary sum g1
real(8) :: x1 ! new g1 value of x
real(8) :: err ! error in convergence of solution
real(8) :: w ! overrelaxation parameter
type(Vector) :: tmpx ! temporary solution vector
! Set overrelaxation parameter
w = ONE
! Dimensions
ng = 1
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
n = A % n
! Perform Gauss Seidel iterations
GS: do igs = 1, 10000
! Check for max iterations met
if (igs == 10000) then
call fatal_error('Maximum Gauss-Seidel iterations encountered.')
endif
! Copy over x vector
call tmpx % copy(x)
! Perform red/black gs iterations
REDBLACK: do irb = 0,1
! Begin loop around matrix rows
ROWS: do irow = 1, n
! Get spatial location
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! Filter out black cells (even)
if (mod(i+j+k,2) == irb) cycle
! Get the index of the diagonals for both rows
call A % search_indices(irow, irow, didx, found)
! Perform temporary sums, first do left of diag block, then right of diag block
tmp1 = ZERO
do icol = A % get_row(irow), didx - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
do icol = didx + 1, A % get_row(irow + 1) - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
! Solve for new x
x1 = (b % val(irow) - tmp1)/A % val(didx)
! Perform overrelaxation
x % val(irow) = (ONE - w)*x % val(irow) + w*x1
end do ROWS
end do REDBLACK
! Check convergence
err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
its = igs
if (err < tol) exit
! Calculation new overrelaxation parameter
w = ONE/(ONE - 0.25_8*cmfd_spectral*w)
end do GS
call tmpx % destroy()
end subroutine cmfd_linsolver_1g
!===============================================================================
! CMFD_LINSOLVER_2G solves the CMFD linear system
!===============================================================================
subroutine cmfd_linsolver_2g(A, b, x, tol, its)
use constants, only: ONE, ZERO
use error, only: fatal_error
use cmfd_header, only: cmfd, cmfd_spectral
type(Matrix), intent(inout) :: A ! coefficient matrix
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! unknown vector
real(8), intent(in) :: tol ! tolerance on final error
integer, intent(out) :: its ! number of inner iterations
integer :: g ! group index
integer :: i ! loop counter for x
integer :: j ! loop counter for y
integer :: k ! loop counter for z
integer :: n ! total size of vector
integer :: nx ! maximum dimension in x direction
integer :: ny ! maximum dimension in y direction
integer :: nz ! maximum dimension in z direction
integer :: ng ! number of energy groups
integer :: d1idx ! index of row "1" diagonal
integer :: d2idx ! index of row "2" diagonal
integer :: igs ! Gauss-Seidel iteration counter
integer :: irb ! Red/Black iteration switch
integer :: irow ! row iteration
integer :: icol ! iteration counter over columns
logical :: found ! did we find col
real(8) :: m11 ! block diagonal component 1,1
real(8) :: m12 ! block diagonal component 1,2
real(8) :: m21 ! block diagonal component 2,1
real(8) :: m22 ! block diagonal component 2,2
real(8) :: dm ! determinant of block diagonal
real(8) :: d11 ! inverse component 1,1
real(8) :: d12 ! inverse component 1,2
real(8) :: d21 ! inverse component 2,1
real(8) :: d22 ! inverse component 2,2
real(8) :: tmp1 ! temporary sum g1
real(8) :: tmp2 ! temporary sum g2
real(8) :: x1 ! new g1 value of x
real(8) :: x2 ! new g2 value of x
real(8) :: err ! error in convergence of solution
real(8) :: w ! overrelaxation parameter
type(Vector) :: tmpx ! temporary solution vector
! Set tolerance and overrelaxation parameter
w = ONE
! Dimensions
ng = 2
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
n = A % n
! Perform Gauss Seidel iterations
GS: do igs = 1, 10000
! Check for max iterations met
if (igs == 10000) then
call fatal_error('Maximum Gauss-Seidel iterations encountered.')
endif
! Copy over x vector
call tmpx % copy(x)
! Perform red/black gs iterations
REDBLACK: do irb = 0,1
! Begin loop around matrix rows
ROWS: do irow = 1, n, 2
! Get spatial location
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! Filter out black cells (even)
if (mod(i+j+k,2) == irb) cycle
! Get the index of the diagonals for both rows
call A % search_indices(irow, irow, d1idx, found)
call A % search_indices(irow + 1, irow + 1, d2idx, found)
! Get block diagonal
m11 = A % val(d1idx) ! group 1 diagonal
m12 = A % val(d1idx + 1) ! group 1 right of diagonal (sorted by col)
m21 = A % val(d2idx - 1) ! group 2 left of diagonal (sorted by col)
m22 = A % val(d2idx) ! group 2 diagonal
! Analytically invert the diagonal
dm = m11*m22 - m12*m21
d11 = m22/dm
d12 = -m12/dm
d21 = -m21/dm
d22 = m11/dm
! Perform temporary sums, first do left of diag block, then right of diag block
tmp1 = ZERO
tmp2 = ZERO
do icol = A % get_row(irow), d1idx - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
do icol = A % get_row(irow + 1), d2idx - 2
tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol))
end do
do icol = d1idx + 2, A % get_row(irow + 1) - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
do icol = d2idx + 1, A % get_row(irow + 2) - 1
tmp2 = tmp2 + A % val(icol)*x % val(A % get_col(icol))
end do
! Adjust with RHS vector
tmp1 = b % val(irow) - tmp1
tmp2 = b % val(irow + 1) - tmp2
! Solve for new x
x1 = d11*tmp1 + d12*tmp2
x2 = d21*tmp1 + d22*tmp2
! Perform overrelaxation
x % val(irow) = (ONE - w)*x % val(irow) + w*x1
x % val(irow + 1) = (ONE - w)*x % val(irow + 1) + w*x2
end do ROWS
end do REDBLACK
! Check convergence
err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
its = igs
if (err < tol) exit
! Calculation new overrelaxation parameter
w = ONE/(ONE - 0.25_8*cmfd_spectral*w)
end do GS
call tmpx % destroy()
end subroutine cmfd_linsolver_2g
!===============================================================================
! CMFD_LINSOLVER_ng solves the CMFD linear system
!===============================================================================
subroutine cmfd_linsolver_ng(A, b, x, tol, its)
use constants, only: ONE, ZERO
use error, only: fatal_error
use cmfd_header, only: cmfd, cmfd_spectral
type(Matrix), intent(inout) :: A ! coefficient matrix
type(Vector), intent(inout) :: b ! right hand side vector
type(Vector), intent(inout) :: x ! unknown vector
real(8), intent(in) :: tol ! tolerance on final error
integer, intent(out) :: its ! number of inner iterations
integer :: g ! group index
integer :: i ! loop counter for x
integer :: j ! loop counter for y
integer :: k ! loop counter for z
integer :: n ! total size of vector
integer :: nx ! maximum dimension in x direction
integer :: ny ! maximum dimension in y direction
integer :: nz ! maximum dimension in z direction
integer :: ng ! number of energy groups
integer :: igs ! Gauss-Seidel iteration counter
integer :: irow ! row iteration
integer :: icol ! iteration counter over columns
integer :: didx ! index for diagonal component
logical :: found ! did we find col
real(8) :: tmp1 ! temporary sum g1
real(8) :: x1 ! new g1 value of x
real(8) :: err ! error in convergence of solution
real(8) :: w ! overrelaxation parameter
type(Vector) :: tmpx ! temporary solution vector
! Set overrelaxation parameter
w = ONE
! Dimensions
ng = 1
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
n = A % n
! Perform Gauss Seidel iterations
GS: do igs = 1, 10000
! Check for max iterations met
if (igs == 10000) then
call fatal_error('Maximum Gauss-Seidel iterations encountered.')
endif
! Copy over x vector
call tmpx % copy(x)
! Begin loop around matrix rows
ROWS: do irow = 1, n
! Get spatial location
call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
! Get the index of the diagonals for both rows
call A % search_indices(irow, irow, didx, found)
! Perform temporary sums, first do left of diag block, then right of diag block
tmp1 = ZERO
do icol = A % get_row(irow), didx - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
do icol = didx + 1, A % get_row(irow + 1) - 1
tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol))
end do
! Solve for new x
x1 = (b % val(irow) - tmp1)/A % val(didx)
! Perform overrelaxation
x % val(irow) = (ONE - w)*x % val(irow) + w*x1
end do ROWS
! Check convergence
err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n)
its = igs
if (err < tol) exit
! Calculation new overrelaxation parameter
w = ONE/(ONE - 0.25_8*cmfd_spectral*w)
end do GS
call tmpx % destroy()
end subroutine cmfd_linsolver_ng
!===============================================================================
! EXTRACT_RESULTS takes results and puts them in CMFD global data object
!===============================================================================
subroutine extract_results()
use cmfd_header, only: cmfd, cmfd_write_matrices
use simulation_header, only: current_batch
character(len=25) :: filename ! name of file to write data
integer :: n ! problem size
! Get problem size
n = loss % n
! Allocate in cmfd object if not already allocated
if (adjoint_calc) then
if (.not. allocated(cmfd%adj_phi)) allocate(cmfd%adj_phi(n))
else
if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n))
end if
! Save values
if (adjoint_calc) then
cmfd % adj_phi = phi_n % val
else
cmfd % phi = phi_n % val
end if
! Save eigenvalue
if(adjoint_calc) then
cmfd%adj_keff = k_n
else
cmfd%keff = k_n
end if
! Normalize phi to 1
if (adjoint_calc) then
cmfd%adj_phi = cmfd%adj_phi/sqrt(sum(cmfd%adj_phi*cmfd%adj_phi))
else
cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi))
end if
! Save dominance ratio
cmfd % dom(current_batch) = norm_n/norm_o
! Write out results
if (cmfd_write_matrices) then
if (adjoint_calc) then
filename = 'adj_fluxvec.dat'
else
filename = 'fluxvec.dat'
end if
call phi_n % write(filename)
end if
end subroutine extract_results
!===============================================================================
! MATRIX_TO_INDICES converts a matrix index to spatial and group indicies
!===============================================================================
subroutine matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz)
use cmfd_header, only: cmfd, cmfd_coremap
integer, intent(out) :: i ! iteration counter for x
integer, intent(out) :: j ! iteration counter for y
integer, intent(out) :: k ! iteration counter for z
integer, intent(out) :: g ! iteration counter for groups
integer, intent(in) :: irow ! iteration counter over row (0 reference)
integer, intent(in) :: nx ! maximum number of x cells
integer, intent(in) :: ny ! maximum number of y cells
integer, intent(in) :: nz ! maximum number of z cells
integer, intent(in) :: ng ! maximum number of groups
! Check for core map
if (cmfd_coremap) then
! Get indices from indexmap
g = mod(irow-1, ng) + 1
i = cmfd % indexmap((irow-1)/ng+1,1)
j = cmfd % indexmap((irow-1)/ng+1,2)
k = cmfd % indexmap((irow-1)/ng+1,3)
else
! Compute indices
g = mod(irow-1, ng) + 1
i = mod(irow-1, ng*nx)/ng + 1
j = mod(irow-1, ng*nx*ny)/(ng*nx)+ 1
k = mod(irow-1, ng*nx*ny*nz)/(ng*nx*ny) + 1
end if
end subroutine matrix_to_indices
!===============================================================================
! FINALIZE frees all memory associated with power iteration
!===============================================================================
subroutine finalize()
! Destroy all objects
call loss % destroy()
call prod % destroy()
call phi_n % destroy()
call phi_o % destroy()
call s_n % destroy()
call s_o % destroy()
call serr_v % destroy
end subroutine finalize
end module cmfd_solver

364
src/cmfd_solver.cpp Normal file
View file

@ -0,0 +1,364 @@
#include <vector>
#include <cmath>
#include "xtensor/xtensor.hpp"
#include "openmc/error.h"
#include "openmc/constants.h"
#include "openmc/capi.h"
namespace openmc {
namespace cmfd {
//==============================================================================
// Global variables
//==============================================================================
std::vector<int> indptr;
std::vector<int> indices;
int dim;
double spectral;
int nx, ny, nz, ng;
xt::xtensor<int, 2> indexmap;
} // namespace cmfd
//==============================================================================
// MATRIX_TO_INDICES converts a matrix index to spatial and group
// indices
//==============================================================================
void matrix_to_indices(int irow, int& g, int& i, int& j, int& k)
{
g = irow % cmfd::ng;
i = cmfd::indexmap(irow/cmfd::ng, 0);
j = cmfd::indexmap(irow/cmfd::ng, 1);
k = cmfd::indexmap(irow/cmfd::ng, 2);
}
//==============================================================================
// GET_DIAGONAL_INDEX returns the index in CSR index array corresponding to
// the diagonal element of a specified row
//==============================================================================
int get_diagonal_index(int row)
{
for (int j = cmfd::indptr[row]; j < cmfd::indptr[row+1]; j++) {
if (cmfd::indices[j] == row)
return j;
}
// Return -1 if not found
return -1;
}
//==============================================================================
// SET_INDEXMAP sets the elements of indexmap based on input coremap
//==============================================================================
void set_indexmap(const int* coremap)
{
for (int z = 0; z < cmfd::nz; z++) {
for (int y = 0; y < cmfd::ny; y++) {
for (int x = 0; x < cmfd::nx; x++) {
if (coremap[(z*cmfd::ny*cmfd::nx) + (y*cmfd::nx) + x] != CMFD_NOACCEL) {
int counter = coremap[(z*cmfd::ny*cmfd::nx) + (y*cmfd::nx) + x];
cmfd::indexmap(counter, 0) = x;
cmfd::indexmap(counter, 1) = y;
cmfd::indexmap(counter, 2) = z;
}
}
}
}
}
//==============================================================================
// CMFD_LINSOLVER_1G solves a one group CMFD linear system
//==============================================================================
int cmfd_linsolver_1g(const double* A_data, const double* b, double* x,
double tol)
{
// Set overrelaxation parameter
double w = 1.0;
// Perform Gauss-Seidel iterations
for (int igs = 1; igs <= 10000; igs++) {
double err = 0.0;
// Copy over x vector
std::vector<double> tmpx {x, x+cmfd::dim};
// Perform red/black Gauss-Seidel iterations
for (int irb = 0; irb < 2; irb++) {
// Loop around matrix rows
for (int irow = 0; irow < cmfd::dim; irow++) {
int g, i, j, k;
matrix_to_indices(irow, g, i, j, k);
// Filter out black cells
if ((i+j+k) % 2 != irb) continue;
// Get index of diagonal for current row
int didx = get_diagonal_index(irow);
// Perform temporary sums, first do left of diag, then right of diag
double tmp1 = 0.0;
for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
// Solve for new x
double x1 = (b[irow] - tmp1) / A_data[didx];
// Perform overrelaxation
x[irow] = (1.0 - w) * x[irow] + w * x1;
// Compute residual and update error
double res = (tmpx[irow] - x[irow]) / tmpx[irow];
err += res * res;
}
}
// Check convergence
err = std::sqrt(err / cmfd::dim);
if (err < tol)
return igs;
// Calculate new overrelaxation parameter
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
}
// Throw error, as max iterations met
fatal_error("Maximum Gauss-Seidel iterations encountered.");
// Return -1 by default, although error thrown before reaching this point
return -1;
}
//==============================================================================
// CMFD_LINSOLVER_2G solves a two group CMFD linear system
//==============================================================================
int cmfd_linsolver_2g(const double* A_data, const double* b, double* x,
double tol)
{
// Set overrelaxation parameter
double w = 1.0;
// Perform Gauss-Seidel iterations
for (int igs = 1; igs <= 10000; igs++) {
double err = 0.0;
// Copy over x vector
std::vector<double> tmpx {x, x+cmfd::dim};
// Perform red/black Gauss-Seidel iterations
for (int irb = 0; irb < 2; irb++) {
// Loop around matrix rows
for (int irow = 0; irow < cmfd::dim; irow+=2) {
int g, i, j, k;
matrix_to_indices(irow, g, i, j, k);
// Filter out black cells
if ((i+j+k) % 2 != irb) continue;
// Get index of diagonals for current row and next row
int d1idx = get_diagonal_index(irow);
int d2idx = get_diagonal_index(irow+1);
// Get block diagonal
double m11 = A_data[d1idx]; // group 1 diagonal
double m12 = A_data[d1idx + 1]; // group 1 right of diagonal (sorted by col)
double m21 = A_data[d2idx - 1]; // group 2 left of diagonal (sorted by col)
double m22 = A_data[d2idx]; // group 2 diagonal
// Analytically invert the diagonal
double dm = m11*m22 - m12*m21;
double d11 = m22/dm;
double d12 = -m12/dm;
double d21 = -m21/dm;
double d22 = m11/dm;
// Perform temporary sums, first do left of diag, then right of diag
double tmp1 = 0.0;
double tmp2 = 0.0;
for (int icol = cmfd::indptr[irow]; icol < d1idx; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
for (int icol = cmfd::indptr[irow+1]; icol < d2idx-1; icol++)
tmp2 += A_data[icol] * x[cmfd::indices[icol]];
for (int icol = d1idx + 2; icol < cmfd::indptr[irow + 1]; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
for (int icol = d2idx + 1; icol < cmfd::indptr[irow + 2]; icol++)
tmp2 += A_data[icol] * x[cmfd::indices[icol]];
// Adjust with RHS vector
tmp1 = b[irow] - tmp1;
tmp2 = b[irow + 1] - tmp2;
// Solve for new x
double x1 = d11*tmp1 + d12*tmp2;
double x2 = d21*tmp1 + d22*tmp2;
// Perform overrelaxation
x[irow] = (1.0 - w) * x[irow] + w * x1;
x[irow + 1] = (1.0 - w) * x[irow + 1] + w * x2;
// Compute residual and update error
double res = (tmpx[irow] - x[irow]) / tmpx[irow];
err += res * res;
}
}
// Check convergence
err = std::sqrt(err / cmfd::dim);
if (err < tol)
return igs;
// Calculate new overrelaxation parameter
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
}
// Throw error, as max iterations met
fatal_error("Maximum Gauss-Seidel iterations encountered.");
// Return -1 by default, although error thrown before reaching this point
return -1;
}
//==============================================================================
// CMFD_LINSOLVER_NG solves a general CMFD linear system
//==============================================================================
int cmfd_linsolver_ng(const double* A_data, const double* b, double* x,
double tol)
{
// Set overrelaxation parameter
double w = 1.0;
// Perform Gauss-Seidel iterations
for (int igs = 1; igs <= 10000; igs++) {
double err = 0.0;
// Copy over x vector
std::vector<double> tmpx {x, x+cmfd::dim};
// Loop around matrix rows
for (int irow = 0; irow < cmfd::dim; irow++) {
// Get index of diagonal for current row
int didx = get_diagonal_index(irow);
// Perform temporary sums, first do left of diag, then right of diag
double tmp1 = 0.0;
for (int icol = cmfd::indptr[irow]; icol < didx; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
for (int icol = didx + 1; icol < cmfd::indptr[irow + 1]; icol++)
tmp1 += A_data[icol] * x[cmfd::indices[icol]];
// Solve for new x
double x1 = (b[irow] - tmp1) / A_data[didx];
// Perform overrelaxation
x[irow] = (1.0 - w) * x[irow] + w * x1;
// Compute residual and update error
double res = (tmpx[irow] - x[irow]) / tmpx[irow];
err += res * res;
}
// Check convergence
err = std::sqrt(err / cmfd::dim);
if (err < tol)
return igs;
// Calculate new overrelaxation parameter
w = 1.0/(1.0 - 0.25 * cmfd::spectral * w);
}
// Throw error, as max iterations met
fatal_error("Maximum Gauss-Seidel iterations encountered.");
// Return -1 by default, although error thrown before reaching this point
return -1;
}
//==============================================================================
// OPENMC_INITIALIZE_LINSOLVER sets the fixed variables that are used for the
// linear solver
//==============================================================================
extern "C"
void openmc_initialize_linsolver(const int* indptr, int len_indptr,
const int* indices, int n_elements, int dim,
double spectral, const int* cmfd_indices,
const int* map)
{
// Store elements of indptr
for (int i = 0; i < len_indptr; i++)
cmfd::indptr.push_back(indptr[i]);
// Store elements of indices
for (int i = 0; i < n_elements; i++)
cmfd::indices.push_back(indices[i]);
// Set dimenion of CMFD problem and specral radius
cmfd::dim = dim;
cmfd::spectral = spectral;
// Set number of groups
cmfd::ng = cmfd_indices[3];
// Set problem dimensions and indexmap if 1 or 2 group problem
if (cmfd::ng == 1 || cmfd::ng == 2) {
cmfd::nx = cmfd_indices[0];
cmfd::ny = cmfd_indices[1];
cmfd::nz = cmfd_indices[2];
// Resize indexmap and set its elements
cmfd::indexmap.resize({static_cast<size_t>(dim), 3});
set_indexmap(map);
}
}
//==============================================================================
// OPENMC_RUN_LINSOLVER runs a Gauss Seidel linear solver to solve CMFD matrix
// equations
//==============================================================================
extern "C"
int openmc_run_linsolver(const double* A_data, const double* b, double* x,
double tol)
{
switch (cmfd::ng) {
case 1:
return cmfd_linsolver_1g(A_data, b, x, tol);
case 2:
return cmfd_linsolver_2g(A_data, b, x, tol);
default:
return cmfd_linsolver_ng(A_data, b, x, tol);
}
}
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" void free_memory_cmfd()
{
cmfd::indptr.clear();
cmfd::indices.clear();
// Resize indexmap to be an empty array
cmfd::indexmap.resize({0});
}
} // namespace openmc

12
src/constants.F90
View file

@ -439,18 +439,6 @@ module constants
ELECTRON_LED = 1, & ! Local Energy Deposition
ELECTRON_TTB = 2 ! Thick Target Bremsstrahlung
!=============================================================================
! CMFD CONSTANTS
! for non-accelerated regions on coarse mesh overlay
integer, parameter :: CMFD_NOACCEL = 99999
! constant to represent a zero flux "albedo"
real(8), parameter :: ZERO_FLUX = 999.0_8
! constant for writing out no residual
real(8), parameter :: CMFD_NORES = 99999.0_8
!=============================================================================
! DELAYED NEUTRON PRECURSOR CONSTANTS

4
src/input_xml.F90
View file

@ -3,8 +3,6 @@ module input_xml
use, intrinsic :: ISO_C_BINDING
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use cmfd_header, only: index_cmfd_mesh
use constants
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use endf, only: reaction_name
@ -153,8 +151,6 @@ contains
! Initialize distribcell_filters
call prepare_distribcell()
if (cmfd_run) call configure_cmfd()
if (run_mode == MODE_PLOTTING) then
! Read plots.xml if it exists
call read_plots_xml()

44
src/output.F90
View file

@ -3,7 +3,6 @@ module output
use, intrinsic :: ISO_C_BINDING
use, intrinsic :: ISO_FORTRAN_ENV
use cmfd_header
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
@ -290,30 +289,12 @@ contains
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " k "
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "Entropy "
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') " Average k "
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') " CMFD k "
select case(trim(cmfd_display))
case('entropy')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "CMFD Ent"
case('balance')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "RMS Bal "
case('source')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "RMS Src "
case('dominance')
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "Dom Rat "
end select
end if
write(UNIT=ou, FMT=*)
write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "========="
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
write(UNIT=ou, FMT='(A20,3X)', ADVANCE='NO') "===================="
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (cmfd_display /= '') &
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
end if
write(UNIT=ou, FMT=*)
end subroutine print_columns
@ -386,26 +367,6 @@ contains
write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO')
end if
! write out cmfd keff if it is active and other display info
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % k_cmfd(current_batch)
select case(trim(cmfd_display))
case('entropy')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % entropy(current_batch)
case('balance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % balance(current_batch)
case('source')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % src_cmp(current_batch)
case('dominance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % dom(current_batch)
end select
end if
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
@ -442,11 +403,6 @@ contains
write(ou,100) " SEND/RECV source sites", time_bank_sendrecv_elapsed()
end if
write(ou,100) " Time accumulating tallies", time_tallies_elapsed()
if (cmfd_run) write(ou,100) " Time in CMFD", time_cmfd % elapsed
if (cmfd_run) write(ou,100) " Building matrices", &
time_cmfdbuild % elapsed
if (cmfd_run) write(ou,100) " Solving matrices", &
time_cmfdsolve % elapsed
write(ou,100) "Total time for finalization", time_finalize_elapsed()
write(ou,100) "Total time elapsed", time_total_elapsed()

53
src/relaxng/cmfd.rnc
View file

@ -1,53 +0,0 @@
element cmfd {
element mesh {
(element dimension { list { xsd:int+ } } |
attribute dimension { list { xsd:int+ } }) &
(element lower_left { list { xsd:double+ } } |
attribute lower_left { list { xsd:double+ } }) &
(
(element upper_right { list { xsd:double+ } } |
attribute upper_right { list { xsd:double+ } }) |
(element width { list { xsd:double+ } } |
attribute width { list { xsd:double+ } })
) &
(element albedo { list { xsd:double+ } } |
attribute albedo { list { xsd:double+ } }) &
(element map { list { xsd:int+ } } |
attribute map { list { xsd:int+ } })? &
(element energy { list { xsd:double+ } } |
attribute energy { list { xsd:double+ } })?
} &
element norm { xsd:double }? &
element feedback { xsd:boolean }? &
element downscatter { xsd:boolean }? &
element dhat_reset { xsd:boolean }? &
element power_monitor { xsd:boolean }? &
element write_matrices { xsd:boolean }? &
element run_adjoint { xsd:boolean }? &
element write_hdf5 { xsd:boolean }? &
element begin { xsd:int }? &
element tally_reset { list { xsd:int+ } }? &
element display { xsd:string }? &
element spectral { xsd:double }? &
element shift { xsd:double }? &
element ktol { xsd:double }? &
element stol { xsd:double }? &
element gauss_seidel_tolerance { list { xsd:double+ } }?
}

215
src/relaxng/cmfd.rng
View file

@ -1,215 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<element name="cmfd" xmlns="http://relaxng.org/ns/structure/1.0" datatypeLibrary="http://www.w3.org/2001/XMLSchema-datatypes">
<interleave>
<element name="mesh">
<interleave>
<choice>
<element name="dimension">
<list>
<oneOrMore>
<data type="int"/>
</oneOrMore>
</list>
</element>
<attribute name="dimension">
<list>
<oneOrMore>
<data type="int"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="lower_left">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<choice>
<element name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="upper_right">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<choice>
<element name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="width">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</choice>
<choice>
<element name="albedo">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="albedo">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
<optional>
<choice>
<element name="map">
<list>
<oneOrMore>
<data type="int"/>
</oneOrMore>
</list>
</element>
<attribute name="map">
<list>
<oneOrMore>
<data type="int"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="energy">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
<attribute name="energy">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
<optional>
<element name="norm">
<data type="double"/>
</element>
</optional>
<optional>
<element name="feedback">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="downscatter">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="dhat_reset">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="power_monitor">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="write_matrices">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="run_adjoint">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="write_hdf5">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="begin">
<data type="int"/>
</element>
</optional>
<optional>
<element name="tally_reset">
<list>
<oneOrMore>
<data type="int"/>
</oneOrMore>
</list>
</element>
</optional>
<optional>
<element name="display">
<data type="string"/>
</element>
</optional>
<optional>
<element name="spectral">
<data type="double"/>
</element>
</optional>
<optional>
<element name="shift">
<data type="double"/>
</element>
</optional>
<optional>
<element name="ktol">
<data type="double"/>
</element>
</optional>
<optional>
<element name="stol">
<data type="double"/>
</element>
</optional>
<optional>
<element name="gauss_seidel_tolerance">
<list>
<oneOrMore>
<data type="double"/>
</oneOrMore>
</list>
</element>
</optional>
</interleave>
</element>

2
src/relaxng/settings.rnc
View file

@ -59,8 +59,6 @@ element settings {
element dagmc { xsd:boolean }? &
element run_cmfd { xsd:boolean }? &
element run_mode { xsd:string }? &
element seed { xsd:positiveInteger }? &

5
src/relaxng/settings.rng
View file

@ -259,11 +259,6 @@
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="run_cmfd">
<data type="boolean"/>
</element>
</optional>
<optional>
<element name="run_mode">
<data type="string"/>

3
src/settings.F90
View file

@ -115,9 +115,6 @@ module settings
real(C_DOUBLE), bind(C) :: res_scat_energy_min
real(C_DOUBLE), bind(C) :: res_scat_energy_max
! Is CMFD active
logical(C_BOOL), bind(C) :: cmfd_run
! No reduction at end of batch
logical(C_BOOL), bind(C) :: reduce_tallies

5
src/settings.cpp
View file

@ -689,11 +689,6 @@ void read_settings_xml()
}
}
// Check for cmfd run
if (check_for_node(root, "run_cmfd")) {
cmfd_run = get_node_value_bool(root, "run_cmfd");
}
// Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) {
xml_node node_res_scat = root.child("resonance_scattering");

15
src/simulation.cpp
View file

@ -25,14 +25,9 @@
namespace openmc {
// data/functions from Fortran side
extern "C" bool cmfd_on;
extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void cmfd_init_batch();
extern "C" void cmfd_tally_init();
extern "C" void execute_cmfd();
extern "C" void init_tally_routines();
extern "C" void load_state_point();
extern "C" void print_batch_keff();
@ -96,9 +91,6 @@ int openmc_simulation_init()
// Allocate tally results arrays if they're not allocated yet
allocate_tally_results();
// Activate the CMFD tallies
cmfd_tally_init();
// Call Fortran initialization
simulation_init_f();
set_micro_xs();
@ -349,11 +341,6 @@ void initialize_batch()
}
}
// check CMFD initialize batch
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
if (settings::cmfd_run) cmfd_init_batch();
}
// Add user tallies to active tallies list
setup_active_tallies();
}
@ -373,8 +360,6 @@ void finalize_batch()
}
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// Perform CMFD calculation if on
if (cmfd_on) execute_cmfd();
// Write batch output
if (mpi::master && settings::verbosity >= 7) print_batch_keff();
}

48
src/state_point.F90
View file

@ -14,7 +14,6 @@ module state_point
use, intrinsic :: ISO_C_BINDING
use bank_header, only: Bank
use cmfd_header
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
@ -62,7 +61,7 @@ contains
integer :: i_xs
integer, allocatable :: id_array(:)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group, tallies_group, tally_group, &
integer(HID_T) :: tallies_group, tally_group, &
filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
character(MAX_WORD_LEN), allocatable :: str_array(:)
@ -167,23 +166,6 @@ contains
! Write out information for eigenvalue run
if (run_mode == MODE_EIGENVALUE) then
call write_eigenvalue_hdf5(file_id)
! Write out CMFD info
if (cmfd_on) then
call write_attribute(file_id, "cmfd_on", 1)
cmfd_group = create_group(file_id, "cmfd")
call write_dataset(cmfd_group, "indices", cmfd % indices)
call write_dataset(cmfd_group, "k_cmfd", cmfd % k_cmfd)
call write_dataset(cmfd_group, "cmfd_src", cmfd % cmfd_src)
call write_dataset(cmfd_group, "cmfd_entropy", cmfd % entropy)
call write_dataset(cmfd_group, "cmfd_balance", cmfd % balance)
call write_dataset(cmfd_group, "cmfd_dominance", cmfd % dom)
call write_dataset(cmfd_group, "cmfd_srccmp", cmfd % src_cmp)
call close_group(cmfd_group)
else
call write_attribute(file_id, "cmfd_on", 0)
end if
end if
tallies_group = create_group(file_id, "tallies")
@ -411,13 +393,6 @@ contains
end if
call write_dataset(runtime_group, "accumulating tallies", &
time_tallies_elapsed())
if (cmfd_run) then
call write_dataset(runtime_group, "CMFD", time_cmfd % get_value())
call write_dataset(runtime_group, "CMFD building matrices", &
time_cmfdbuild % get_value())
call write_dataset(runtime_group, "CMFD solving matrices", &
time_cmfdsolve % get_value())
end if
call write_dataset(runtime_group, "total", time_total_elapsed())
call close_group(runtime_group)
@ -451,7 +426,6 @@ contains
integer, allocatable :: array(:)
integer(C_INT64_T) :: seed
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
logical :: source_present
@ -546,26 +520,6 @@ contains
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))
! Read in to see if CMFD was on
call read_attribute(int_array(1), file_id, "cmfd_on")
! Read in CMFD info
if (int_array(1) == 1) then
cmfd_group = open_group(file_id, "cmfd")
call read_dataset(cmfd % indices, cmfd_group, "indices")
call read_dataset(cmfd % k_cmfd(1:restart_batch), cmfd_group, "k_cmfd")
call read_dataset(cmfd % cmfd_src, cmfd_group, "cmfd_src")
call read_dataset(cmfd % entropy(1:restart_batch), cmfd_group, &
"cmfd_entropy")
call read_dataset(cmfd % balance(1:restart_batch), cmfd_group, &
"cmfd_balance")
call read_dataset(cmfd % dom(1:restart_batch), cmfd_group, &
"cmfd_dominance")
call read_dataset(cmfd % src_cmp(1:restart_batch), cmfd_group, &
"cmfd_srccmp")
call close_group(cmfd_group)
end if
end if
! Read number of realizations for global tallies

15
tests/regression_tests/cmfd_feed/cmfd.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<cmfd>
<mesh>
<lower_left> -10 -1 -1 </lower_left>
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
</mesh>
<begin>5</begin>
<display>dominance</display>
<feedback>true</feedback>
<gauss_seidel_tolerance>1.e-15 1.e-20</gauss_seidel_tolerance>
</cmfd>

34
tests/regression_tests/cmfd_feed/results_true.dat
View file

@ -391,10 +391,6 @@ cmfd indices
1.000000E+00
1.000000E+00
k cmfd
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.170416E+00
1.172966E+00
1.165537E+00
@ -412,10 +408,6 @@ k cmfd
1.167848E+00
1.165116E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.203643E+00
3.207943E+00
3.213367E+00
@ -433,10 +425,6 @@ cmfd entropy
3.233193E+00
3.232564E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.009062E-03
4.431773E-03
3.152666E-03
@ -446,18 +434,14 @@ cmfd balance
1.502427E-03
1.589825E-03
1.566020E-03
1.219160E-03
1.017888E-03
9.771622E-04
1.010120E-03
1.073382E-03
1.172758E-03
9.827332E-04
1.170009E-03
9.507594E-04
9.072590E-04
1.009003E-03
1.064703E-03
1.163610E-03
9.756308E-04
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
5.397E-01
5.425E-01
5.481E-01
@ -475,10 +459,6 @@ cmfd dominance ratio
5.524E-01
5.523E-01
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.959834E-03
5.655657E-03
3.886185E-03

3
tests/regression_tests/cmfd_feed/settings.xml
View file

@ -23,7 +23,4 @@
</mesh>
<entropy_mesh>10</entropy_mesh>
<!-- Run CMFD -->
<run_cmfd>true</run_cmfd>
</settings>

134
tests/regression_tests/cmfd_feed/test.py
View file

@ -1,6 +1,138 @@
from tests.testing_harness import CMFDTestHarness
from openmc import cmfd
import numpy as np
import scipy.sparse
def test_cmfd_physical_adjoint():
"""Test physical adjoint functionality of CMFD
This test runs CMFD with a physical adjoint calculation and asserts that
the adjoint k-effective and flux vector are equal to the non-adjoint
k-effective and flux vector at the last batch (equivalent for 1 group
problems).
"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -10.0, -1.0, -1.0
cmfd_mesh.upper_right = 10.0, 1.0, 1.0
cmfd_mesh.dimension = 10, 1, 1
cmfd_mesh.albedo = 0.0, 0.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run_adjoint = True
cmfd_run.adjoint_type = 'physical'
cmfd_run.run()
assert(np.all(cmfd_run._phi == cmfd_run._adj_phi))
assert(cmfd_run._adj_keff == cmfd_run._keff)
def test_cmfd_math_adjoint():
"""Test mathematical adjoint functionality of CMFD
This test runs CMFD with a mathematical adjoint calculation and asserts
that the adjoint k-effective and flux vector are equal to the non-adjoint
k-effective and flux vector at the last batch (equivalent for 1 group
problems).
"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -10.0, -1.0, -1.0
cmfd_mesh.upper_right = 10.0, 1.0, 1.0
cmfd_mesh.dimension = 10, 1, 1
cmfd_mesh.albedo = 0.0, 0.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run_adjoint = True
cmfd_run.adjoint_type = 'math'
cmfd_run.run()
assert(np.all(cmfd_run._phi == cmfd_run._adj_phi))
assert(cmfd_run._adj_keff == cmfd_run._keff)
def test_cmfd_write_matrices():
"""Test write matrices functionality of CMFD
This test runs CMFD with feedback and loads the loss/production matrices
and flux vector that are saved to disk, and checks to make sure these
values are consistent with each other and simulation results.
"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -10.0, -1.0, -1.0
cmfd_mesh.upper_right = 10.0, 1.0, 1.0
cmfd_mesh.dimension = 10, 1, 1
cmfd_mesh.albedo = 0.0, 0.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.write_matrices = True
cmfd_run.run()
# Load loss matrix from numpy output file
loss_np = scipy.sparse.load_npz('loss.npz').todense()
# Load loss matrix from data file
loss_dat = np.loadtxt("loss.dat", delimiter=',')
# Go through each element of loss_dat and compare to loss_np
for elem in loss_dat:
assert(np.isclose(loss_np[int(elem[0]), int(elem[1])], elem[2]))
# Load production matrix from numpy output file
prod_np = scipy.sparse.load_npz('prod.npz').todense()
# Load production matrix from data file
prod_dat = np.loadtxt("prod.dat", delimiter=',')
# Go through each element of prod_dat and compare to prod_np
for elem in prod_dat:
assert(np.isclose(prod_np[int(elem[0]), int(elem[1])], elem[2]))
# Load flux vector from numpy output file
flux_np = np.load('fluxvec.npy')
# Load flux from data file
flux_dat = np.loadtxt("fluxvec.dat", delimiter='\n')
# Compare flux from numpy file, .dat file, and from simulation
assert(np.all(np.isclose(flux_np, cmfd_run._phi)))
assert(np.all(np.isclose(flux_np, flux_dat)))
def test_cmfd_feed():
harness = CMFDTestHarness('statepoint.20.h5')
"""Test 1 group CMFD solver with CMFD feedback"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -10.0, -1.0, -1.0
cmfd_mesh.upper_right = 10.0, 1.0, 1.0
cmfd_mesh.dimension = 10, 1, 1
cmfd_mesh.albedo = 0.0, 0.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run()
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()

0
tests/regression_tests/cmfd_feed_2g/__init__.py Normal file
View file

31
tests/regression_tests/cmfd_feed_2g/geometry.xml Normal file
View file

@ -0,0 +1,31 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" fill="5" region="151 -152 153 -154 1380 -1370" />
<cell id="52" material="9" name="Fuel rod active region - 1.6% enr radial 0: Fuel 1.6%" region="-35" universe="8" />
<cell id="53" material="5" name="Fuel rod active region - 1.6% enr radial 1: Helium" region="35 -36" universe="8" />
<cell id="54" material="7" name="Fuel rod active region - 1.6% enr radial outer: Zircaloy 4" region="36" universe="8" />
<cell fill="8" id="1253" name="Fuel cell" region="-150" universe="763" />
<cell id="1254" material="15" name="Water" region="150" universe="763" />
<lattice id="5">
<dimension>2 2</dimension>
<lower_left>-1.25984 -1.25984</lower_left>
<pitch>1.25984 1.25984</pitch>
<universes>
763 763
763 763
</universes>
</lattice>
<surface coeffs="0.0 0.0 0.39218" id="35" name="Fuel pellet OR" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.40005" id="36" name="Fuel clad IR" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.4572" id="150" name="Fuel clad OR" type="z-cylinder" />
<surface boundary="reflective" coeffs="-1.25984" id="151" name="minimum x" type="x-plane" />
<surface boundary="reflective" coeffs="1.25984" id="152" name="maximum x" type="x-plane" />
<surface boundary="reflective" coeffs="-1.25984" id="153" name="minimum y" type="y-plane" />
<surface boundary="reflective" coeffs="1.25984" id="154" name="maximum y" type="y-plane" />
<surface boundary="reflective" coeffs="1.0" id="1370" name="Highest Extent" type="z-plane" />
<surface boundary="reflective" coeffs="-1.0" id="1380" name="Lowest Extent" type="z-plane" />
</geometry>

59
tests/regression_tests/cmfd_feed_2g/materials.xml Normal file
View file

@ -0,0 +1,59 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="5" name="Helium">
<temperature>300</temperature>
<density units="g/cc" value="0.0015981" />
<nuclide name="He3" wo="1.5070346049256974e-06" />
<nuclide name="He4" wo="0.999998492965395" />
</material>
<material id="7" name="Zircaloy 4">
<temperature>300</temperature>
<density units="g/cc" value="6.55" />
<nuclide name="O16" wo="0.0012494965182849112" />
<nuclide name="O17" wo="5.034817150887735e-07" />
<nuclide name="Cr50" wo="4.1736864731106146e-05" />
<nuclide name="Cr52" wo="0.0008369936242576807" />
<nuclide name="Cr53" wo="9.673586881507429e-05" />
<nuclide name="Cr54" wo="2.4533642196138756e-05" />
<nuclide name="Fe54" wo="0.00011855672274761877" />
<nuclide name="Fe56" wo="0.001929932104229657" />
<nuclide name="Fe57" wo="4.536774095388075e-05" />
<nuclide name="Fe58" wo="6.143432068843669e-06" />
<nuclide name="Zr90" wo="0.49750307249921255" />
<nuclide name="Zr91" wo="0.10970127796055709" />
<nuclide name="Zr92" wo="0.16952409354767467" />
<nuclide name="Zr94" wo="0.17553856942304608" />
<nuclide name="Zr96" wo="0.02888298656950975" />
<nuclide name="Sn112" wo="0.0001325869644430062" />
<nuclide name="Sn114" wo="9.182449637587617e-05" />
<nuclide name="Sn115" wo="4.771905922545867e-05" />
<nuclide name="Sn116" wo="0.002058423153629443" />
<nuclide name="Sn117" wo="0.0010966473429083066" />
<nuclide name="Sn118" wo="0.0034879812938438245" />
<nuclide name="Sn119" wo="0.001247577110245757" />
<nuclide name="Sn120" wo="0.004771539495238715" />
<nuclide name="Sn122" wo="0.0006894094798456136" />
<nuclide name="Sn124" wo="0.000876291604244001" />
</material>
<material depletable="true" id="9" name="Fuel 1.6%">
<temperature>300</temperature>
<density units="g/cc" value="10.31341" />
<nuclide ao="1.9992419999999993" name="O16" />
<nuclide ao="0.0007579999999999998" name="O17" />
<nuclide ao="0.00014571996318473284" name="U234" />
<nuclide ao="0.01630316269333107" name="U235" />
<nuclide ao="0.983476441027697" name="U238" />
<nuclide ao="7.467631578722921e-05" name="U236" />
</material>
<material id="15" name="Borated Water">
<temperature>300</temperature>
<density units="g/cc" value="0.7405820675158279" />
<nuclide ao="0.00032178659941803253" name="B10" />
<nuclide ao="0.0013017583017829388" name="B11" />
<nuclide ao="1.996441935899364" name="H1" />
<nuclide ao="0.0003109742982341739" name="H2" />
<nuclide ao="0.9979980704223166" name="O16" />
<nuclide ao="0.0003783846764824448" name="O17" />
<sab name="c_H_in_H2O" />
</material>
</materials>

434
tests/regression_tests/cmfd_feed_2g/results_true.dat Normal file
View file

@ -0,0 +1,434 @@
k-combined:
1.038883E+00 1.017026E-02
tally 1:
1.167304E+02
1.362680E+03
1.157179E+02
1.339273E+03
1.167087E+02
1.362782E+03
1.164493E+02
1.356411E+03
tally 2:
4.400233E+01
9.698989E+01
6.541957E+01
2.144299E+02
1.863245E+02
1.738212E+03
1.038047E+02
5.390768E+02
4.368854E+01
9.564001E+01
6.489620E+01
2.109923E+02
1.836691E+02
1.687792E+03
1.022706E+02
5.232407E+02
4.433402E+01
9.859049E+01
6.548060E+01
2.150470E+02
1.839219E+02
1.692209E+03
1.031634E+02
5.323852E+02
4.380686E+01
9.606723E+01
6.480355E+01
2.103905E+02
1.868847E+02
1.747254E+03
1.038109E+02
5.391360E+02
tally 3:
6.192979E+01
1.921761E+02
0.000000E+00
0.000000E+00
2.033842E-02
4.375294E-05
4.296338E+00
9.280716E-01
3.493776E+00
6.157094E-01
0.000000E+00
0.000000E+00
9.876716E+01
4.880113E+02
9.043140E-01
4.156054E-02
6.138292E+01
1.887853E+02
0.000000E+00
0.000000E+00
1.591722E-02
2.328432E-05
4.296798E+00
9.289349E-01
3.584967E+00
6.444749E-01
0.000000E+00
0.000000E+00
9.717439E+01
4.724254E+02
8.435691E-01
3.666151E-02
6.192881E+01
1.924016E+02
0.000000E+00
0.000000E+00
1.796382E-02
3.190854E-05
4.343238E+00
9.514039E-01
3.504683E+00
6.157615E-01
0.000000E+00
0.000000E+00
9.818871E+01
4.822799E+02
8.401346E-01
3.664037E-02
6.130607E+01
1.882837E+02
0.000000E+00
0.000000E+00
1.220252E-02
1.711667E-05
4.245413E+00
9.056738E-01
3.468894E+00
6.033679E-01
0.000000E+00
0.000000E+00
9.885292E+01
4.888720E+02
9.509199E-01
4.670951E-02
tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.200260E+00
4.243011E+00
3.728022E+01
6.953801E+01
9.178009E+00
4.218903E+00
3.738925E+01
6.993833E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.019225E+00
4.079528E+00
3.709185E+01
6.883085E+01
9.037478E+00
4.095846E+00
3.710303E+01
6.888447E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.178009E+00
4.218903E+00
3.738925E+01
6.993833E+01
9.200260E+00
4.243011E+00
3.728022E+01
6.953801E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.097360E+00
4.147125E+00
3.700603E+01
6.851110E+01
9.003417E+00
4.059372E+00
3.717266E+01
6.912709E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.042233E+00
4.097018E+00
3.711597E+01
6.891834E+01
9.107597E+00
4.161211E+00
3.715168E+01
6.904544E+01
9.037478E+00
4.095846E+00
3.710303E+01
6.888447E+01
9.019225E+00
4.079528E+00
3.709185E+01
6.883085E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.107597E+00
4.161211E+00
3.715168E+01
6.904544E+01
9.042233E+00
4.097018E+00
3.711597E+01
6.891834E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.003417E+00
4.059372E+00
3.717266E+01
6.912709E+01
9.097360E+00
4.147125E+00
3.700603E+01
6.851110E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
tally 5:
6.195013E+01
1.923032E+02
1.022516E+02
5.230716E+02
1.418310E+01
1.011490E+01
4.677537E+01
1.094859E+02
6.139884E+01
1.888834E+02
1.007512E+02
5.078162E+02
1.403236E+01
9.884347E+00
4.598486E+01
1.058070E+02
6.194677E+01
1.925125E+02
1.016864E+02
5.172531E+02
1.407797E+01
9.971785E+00
4.623118E+01
1.069663E+02
6.131828E+01
1.883600E+02
1.023107E+02
5.236668E+02
1.375796E+01
9.501501E+00
4.638162E+01
1.076532E+02
cmfd indices
2.000000E+00
2.000000E+00
1.000000E+00
2.000000E+00
k cmfd
1.034617E+00
1.034521E+00
1.036479E+00
1.035658E+00
1.032505E+00
1.027277E+00
1.029326E+00
1.029881E+00
1.035579E+00
1.034753E+00
1.031963E+00
1.034797E+00
1.035383E+00
1.035908E+00
1.035629E+00
1.036034E+00
cmfd entropy
1.999493E+00
1.999058E+00
1.999331E+00
1.999370E+00
1.999383E+00
1.999455E+00
1.999640E+00
1.999810E+00
1.999881E+00
1.999937E+00
1.999973E+00
1.999965E+00
1.999972E+00
1.999982E+00
1.999998E+00
1.999984E+00
cmfd balance
3.993043E-04
4.168558E-04
6.384693E-04
3.928217E-04
3.789833E-04
2.684860E-04
4.849910E-04
1.084019E-03
1.091775E-03
5.459767E-04
4.455545E-04
4.011468E-04
3.710248E-04
3.577152E-04
3.845993E-04
3.900672E-04
cmfd dominance ratio
6.239E-03
6.303E-03
6.347E-03
6.388E-03
6.340E-03
6.334E-03
6.270E-03
6.173E-03
6.068E-03
6.061E-03
6.016E-03
6.108E-03
6.117E-03
6.074E-03
6.005E-03
5.996E-03
cmfd openmc source comparison
1.931386E-05
2.839162E-05
1.407962E-05
6.718148E-06
9.164193E-06
1.382539E-05
2.871606E-05
4.192300E-05
5.327516E-05
6.566500E-05
6.655541E-05
6.034977E-05
6.004677E-05
5.711623E-05
6.271264E-05
6.425363E-05
cmfd source
2.510278E-01
2.480592E-01
2.501750E-01
2.507380E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

25
tests/regression_tests/cmfd_feed_2g/settings.xml Normal file
View file

@ -0,0 +1,25 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- Starting source -->
<source strength="1.0">
<space type="fission">
<parameters>-1.25984 -1.25984 -1 1.25984 1.25984 1</parameters>
</space>
</source>
<!-- Shannon Entropy -->
<mesh id="10">
<dimension> 2 2 1 </dimension>
<lower_left> -1.25984 -1.25984 -1.0 </lower_left>
<upper_right> 1.25984 1.25984 1.0 </upper_right>
</mesh>
<entropy_mesh>10</entropy_mesh>
</settings>

21
tests/regression_tests/cmfd_feed_2g/tallies.xml Normal file
View file

@ -0,0 +1,21 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>regular</type>
<lower_left>-1.25984 -1.25984 -1.0 </lower_left>
<upper_right>1.25984 1.25984 1.0</upper_right>
<dimension>2 2 1</dimension>
</mesh>
<filter id="1">
<type>mesh</type>
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>

28
tests/regression_tests/cmfd_feed_2g/test.py Normal file
View file

@ -0,0 +1,28 @@
from tests.testing_harness import CMFDTestHarness
from openmc import cmfd
import numpy as np
def test_cmfd_feed_2g():
"""Test 2 group CMFD solver results with CMFD feedback"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -1.25984, -1.25984, -1.0
cmfd_mesh.upper_right = 1.25984, 1.25984, 1.0
cmfd_mesh.dimension = 2, 2, 1
cmfd_mesh.energy = [0.0, 0.625, 20000000]
cmfd_mesh.albedo = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.downscatter = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run()
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()

0
tests/regression_tests/cmfd_feed_ng/__init__.py Normal file
View file

31
tests/regression_tests/cmfd_feed_ng/geometry.xml Normal file
View file

@ -0,0 +1,31 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" fill="5" region="151 -152 153 -154 1380 -1370" />
<cell id="52" material="9" name="Fuel rod active region - 1.6% enr radial 0: Fuel 1.6%" region="-35" universe="8" />
<cell id="53" material="5" name="Fuel rod active region - 1.6% enr radial 1: Helium" region="35 -36" universe="8" />
<cell id="54" material="7" name="Fuel rod active region - 1.6% enr radial outer: Zircaloy 4" region="36" universe="8" />
<cell fill="8" id="1253" name="Fuel cell" region="-150" universe="763" />
<cell id="1254" material="15" name="Water" region="150" universe="763" />
<lattice id="5">
<dimension>2 2</dimension>
<lower_left>-1.25984 -1.25984</lower_left>
<pitch>1.25984 1.25984</pitch>
<universes>
763 763
763 763
</universes>
</lattice>
<surface coeffs="0.0 0.0 0.39218" id="35" name="Fuel pellet OR" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.40005" id="36" name="Fuel clad IR" type="z-cylinder" />
<surface coeffs="0.0 0.0 0.4572" id="150" name="Fuel clad OR" type="z-cylinder" />
<surface boundary="reflective" coeffs="-1.25984" id="151" name="minimum x" type="x-plane" />
<surface boundary="reflective" coeffs="1.25984" id="152" name="maximum x" type="x-plane" />
<surface boundary="reflective" coeffs="-1.25984" id="153" name="minimum y" type="y-plane" />
<surface boundary="reflective" coeffs="1.25984" id="154" name="maximum y" type="y-plane" />
<surface boundary="reflective" coeffs="1.0" id="1370" name="Highest Extent" type="z-plane" />
<surface boundary="reflective" coeffs="-1.0" id="1380" name="Lowest Extent" type="z-plane" />
</geometry>

59
tests/regression_tests/cmfd_feed_ng/materials.xml Normal file
View file

@ -0,0 +1,59 @@
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="5" name="Helium">
<temperature>300</temperature>
<density units="g/cc" value="0.0015981" />
<nuclide name="He3" wo="1.5070346049256974e-06" />
<nuclide name="He4" wo="0.999998492965395" />
</material>
<material id="7" name="Zircaloy 4">
<temperature>300</temperature>
<density units="g/cc" value="6.55" />
<nuclide name="O16" wo="0.0012494965182849112" />
<nuclide name="O17" wo="5.034817150887735e-07" />
<nuclide name="Cr50" wo="4.1736864731106146e-05" />
<nuclide name="Cr52" wo="0.0008369936242576807" />
<nuclide name="Cr53" wo="9.673586881507429e-05" />
<nuclide name="Cr54" wo="2.4533642196138756e-05" />
<nuclide name="Fe54" wo="0.00011855672274761877" />
<nuclide name="Fe56" wo="0.001929932104229657" />
<nuclide name="Fe57" wo="4.536774095388075e-05" />
<nuclide name="Fe58" wo="6.143432068843669e-06" />
<nuclide name="Zr90" wo="0.49750307249921255" />
<nuclide name="Zr91" wo="0.10970127796055709" />
<nuclide name="Zr92" wo="0.16952409354767467" />
<nuclide name="Zr94" wo="0.17553856942304608" />
<nuclide name="Zr96" wo="0.02888298656950975" />
<nuclide name="Sn112" wo="0.0001325869644430062" />
<nuclide name="Sn114" wo="9.182449637587617e-05" />
<nuclide name="Sn115" wo="4.771905922545867e-05" />
<nuclide name="Sn116" wo="0.002058423153629443" />
<nuclide name="Sn117" wo="0.0010966473429083066" />
<nuclide name="Sn118" wo="0.0034879812938438245" />
<nuclide name="Sn119" wo="0.001247577110245757" />
<nuclide name="Sn120" wo="0.004771539495238715" />
<nuclide name="Sn122" wo="0.0006894094798456136" />
<nuclide name="Sn124" wo="0.000876291604244001" />
</material>
<material depletable="true" id="9" name="Fuel 1.6%">
<temperature>300</temperature>
<density units="g/cc" value="10.31341" />
<nuclide ao="1.9992419999999993" name="O16" />
<nuclide ao="0.0007579999999999998" name="O17" />
<nuclide ao="0.00014571996318473284" name="U234" />
<nuclide ao="0.01630316269333107" name="U235" />
<nuclide ao="0.983476441027697" name="U238" />
<nuclide ao="7.467631578722921e-05" name="U236" />
</material>
<material id="15" name="Borated Water">
<temperature>300</temperature>
<density units="g/cc" value="0.7405820675158279" />
<nuclide ao="0.00032178659941803253" name="B10" />
<nuclide ao="0.0013017583017829388" name="B11" />
<nuclide ao="1.996441935899364" name="H1" />
<nuclide ao="0.0003109742982341739" name="H2" />
<nuclide ao="0.9979980704223166" name="O16" />
<nuclide ao="0.0003783846764824448" name="O17" />
<sab name="c_H_in_H2O" />
</material>
</materials>

621
tests/regression_tests/cmfd_feed_ng/results_true.dat Normal file
View file

@ -0,0 +1,621 @@
k-combined:
1.030507E+00 9.667655E-03
tally 1:
1.159445E+02
1.344573E+03
1.159285E+02
1.344429E+03
1.161340E+02
1.349209E+03
1.160434E+02
1.347427E+03
tally 2:
3.482032E+01
7.610041E+01
5.123532E+01
1.648226E+02
1.063498E+01
7.110316E+00
9.001212E+00
5.085589E+00
1.435947E+02
1.322778E+03
7.400822E+01
3.425876E+02
3.415144E+01
7.318850E+01
5.046314E+01
1.597317E+02
1.028893E+01
6.661534E+00
8.653997E+00
4.698660E+00
1.442064E+02
1.338979E+03
7.405073E+01
3.428774E+02
3.431740E+01
7.390622E+01
5.109288E+01
1.639090E+02
1.034894E+01
6.746021E+00
8.657005E+00
4.712481E+00
1.383790E+02
1.200016E+03
7.320761E+01
3.351661E+02
3.426169E+01
7.373156E+01
5.042174E+01
1.596859E+02
1.041331E+01
6.828453E+00
8.729287E+00
4.779221E+00
1.374766E+02
1.182015E+03
7.390660E+01
3.416196E+02
tally 3:
4.842028E+01
1.472522E+02
0.000000E+00
0.000000E+00
1.766768E-02
3.669087E-05
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.489090E+00
7.670426E-01
2.522457E+00
3.988419E-01
0.000000E+00
0.000000E+00
6.323118E+00
2.515554E+00
0.000000E+00
0.000000E+00
1.000000E-03
1.000000E-06
9.772798E-02
7.833336E-04
3.020916E-01
5.962947E-03
0.000000E+00
0.000000E+00
2.669802E+00
4.494602E-01
0.000000E+00
0.000000E+00
6.995726E+01
3.061105E+02
6.245439E-01
2.516241E-02
4.767259E+01
1.425559E+02
0.000000E+00
0.000000E+00
1.511055E-02
2.706989E-05
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.359388E+00
7.162206E-01
2.420435E+00
3.693966E-01
0.000000E+00
0.000000E+00
6.097921E+00
2.335376E+00
0.000000E+00
0.000000E+00
2.973834E-03
2.952905E-06
9.585516E-02
7.950754E-04
3.306195E-01
7.226156E-03
0.000000E+00
0.000000E+00
2.606740E+00
4.273433E-01
0.000000E+00
0.000000E+00
7.004531E+01
3.068103E+02
5.609353E-01
2.054700E-02
4.844434E+01
1.474236E+02
0.000000E+00
0.000000E+00
6.909228E-03
6.822381E-06
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.337728E+00
7.036827E-01
2.418767E+00
3.680337E-01
0.000000E+00
0.000000E+00
6.114057E+00
2.352184E+00
0.000000E+00
0.000000E+00
1.000000E-03
1.000000E-06
1.057559E-01
8.659764E-04
3.268970E-01
6.949165E-03
0.000000E+00
0.000000E+00
2.538042E+00
4.047397E-01
0.000000E+00
0.000000E+00
6.922289E+01
2.996756E+02
6.190458E-01
2.456695E-02
4.757150E+01
1.421640E+02
0.000000E+00
0.000000E+00
8.765775E-03
1.054666E-05
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.436793E+00
7.430145E-01
2.473841E+00
3.843375E-01
0.000000E+00
0.000000E+00
6.104077E+00
2.339598E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.061129E-01
8.936114E-04
3.092546E-01
6.384841E-03
0.000000E+00
0.000000E+00
2.539292E+00
4.045561E-01
0.000000E+00
0.000000E+00
6.994273E+01
3.059846E+02
5.887450E-01
2.235244E-02
tally 4:
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.109363E+00
3.174618E+00
2.099346E+00
2.784777E-01
2.785196E+01
4.850659E+01
7.069021E+00
3.132588E+00
2.149589E+00
2.900131E-01
2.768129E+01
4.790811E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.175941E+00
3.224156E+00
2.135542E+00
2.863424E-01
2.728334E+01
4.655111E+01
7.206774E+00
3.256252E+00
2.075924E+00
2.712640E-01
2.747967E+01
4.722188E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.069021E+00
3.132588E+00
2.149589E+00
2.900131E-01
2.768129E+01
4.790811E+01
7.109363E+00
3.174618E+00
2.099346E+00
2.784777E-01
2.785196E+01
4.850659E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.108601E+00
3.175833E+00
2.103259E+00
2.783890E-01
2.784173E+01
4.846581E+01
7.107753E+00
3.173064E+00
2.087728E+00
2.737368E-01
2.779197E+01
4.828855E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.148041E+00
3.216417E+00
2.124230E+00
2.831606E-01
2.743879E+01
4.708424E+01
7.081749E+00
3.151711E+00
2.062609E+00
2.682691E-01
2.761154E+01
4.768316E+01
7.206774E+00
3.256252E+00
2.075924E+00
2.712640E-01
2.747967E+01
4.722188E+01
7.175941E+00
3.224156E+00
2.135542E+00
2.863424E-01
2.728334E+01
4.655111E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.081749E+00
3.151711E+00
2.062609E+00
2.682691E-01
2.761154E+01
4.768316E+01
7.148041E+00
3.216417E+00
2.124230E+00
2.831606E-01
2.743879E+01
4.708424E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.107753E+00
3.173064E+00
2.087728E+00
2.737368E-01
2.779197E+01
4.828855E+01
7.108601E+00
3.175833E+00
2.103259E+00
2.783890E-01
2.784173E+01
4.846581E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
tally 5:
4.843795E+01
1.473599E+02
8.846575E+00
4.912671E+00
7.292412E+01
3.326224E+02
1.096673E+01
7.549697E+00
4.119614E+00
1.064248E+00
3.323487E+01
6.913295E+01
4.768770E+01
1.426471E+02
8.521329E+00
4.554207E+00
7.297983E+01
3.330439E+02
1.083867E+01
7.393362E+00
3.955093E+00
9.810363E-01
3.309790E+01
6.852254E+01
4.845125E+01
1.474656E+02
8.533824E+00
4.579629E+00
7.208076E+01
3.249387E+02
1.100395E+01
7.627044E+00
3.908798E+00
9.610162E-01
3.287630E+01
6.763801E+01
4.758027E+01
1.422161E+02
8.577918E+00
4.614174E+00
7.278634E+01
3.313433E+02
1.072190E+01
7.216342E+00
3.922303E+00
9.640585E-01
3.334768E+01
6.955922E+01
cmfd indices
2.000000E+00
2.000000E+00
1.000000E+00
3.000000E+00
k cmfd
1.024488E+00
1.022339E+00
1.024720E+00
1.025459E+00
1.029180E+00
1.022688E+00
1.019305E+00
1.020356E+00
1.023384E+00
1.023519E+00
1.026360E+00
cmfd entropy
1.999710E+00
1.999709E+00
1.999704E+00
1.999747E+00
1.999885E+00
1.999911E+00
1.999905E+00
1.999936E+00
1.999873E+00
1.999826E+00
1.999768E+00
cmfd balance
3.179889E-04
4.166064E-04
2.746180E-04
3.345930E-04
5.129605E-04
4.008511E-04
3.628451E-04
3.401727E-04
3.166255E-04
2.642675E-04
2.482640E-04
cmfd dominance ratio
3.737E-03
3.786E-03
3.798E-03
3.789E-03
3.741E-03
3.772E-03
3.793E-03
3.684E-03
3.740E-03
3.764E-03
3.798E-03
cmfd openmc source comparison
3.559244E-05
4.429506E-05
2.739679E-05
2.715448E-05
5.105340E-05
4.069206E-05
4.033877E-05
3.669734E-05
3.540373E-05
2.907247E-05
2.734478E-05
cmfd source
2.564687E-01
2.445733E-01
2.473769E-01
2.515811E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

25
tests/regression_tests/cmfd_feed_ng/settings.xml Normal file
View file

@ -0,0 +1,25 @@
<?xml version="1.0" encoding="UTF-8"?>
<settings>
<!-- Parameters for criticality calculation -->
<run_mode>eigenvalue</run_mode>
<batches>20</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- Starting source -->
<source strength="1.0">
<space type="fission">
<parameters>-1.25984 -1.25984 -1 1.25984 1.25984 1</parameters>
</space>
</source>
<!-- Shannon Entropy -->
<mesh id="10">
<dimension> 2 2 1 </dimension>
<lower_left> -1.25984 -1.25984 -1.0 </lower_left>
<upper_right> 1.25984 1.25984 1.0 </upper_right>
</mesh>
<entropy_mesh>10</entropy_mesh>
</settings>

21
tests/regression_tests/cmfd_feed_ng/tallies.xml Normal file
View file

@ -0,0 +1,21 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1">
<type>regular</type>
<lower_left>-1.25984 -1.25984 -1.0 </lower_left>
<upper_right>1.25984 1.25984 1.0</upper_right>
<dimension>2 2 1</dimension>
</mesh>
<filter id="1">
<type>mesh</type>
<bins>1</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>flux</scores>
</tally>
</tallies>

29
tests/regression_tests/cmfd_feed_ng/test.py Normal file
View file

@ -0,0 +1,29 @@
from tests.testing_harness import CMFDTestHarness
from openmc import cmfd
import numpy as np
def test_cmfd_feed_ng():
"""Test n group CMFD solver with CMFD feedback"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -1.25984, -1.25984, -1.0
cmfd_mesh.upper_right = 1.25984, 1.25984, 1.0
cmfd_mesh.dimension = 2, 2, 1
cmfd_mesh.energy = [0.0, 0.625, 5.53080, 20000000]
cmfd_mesh.albedo = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.reset = [5]
cmfd_run.begin = 10
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = True
cmfd_run.downscatter = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run()
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()

16
tests/regression_tests/cmfd_nofeed/cmfd.xml
View file

@ -1,16 +0,0 @@
<?xml version="1.0"?>
<cmfd>
<mesh>
<lower_left> -10 -1 -1 </lower_left>
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
</mesh>
<begin> 5 </begin>
<display> dominance </display>
<feedback> false </feedback>
<gauss_seidel_tolerance> 1.e-15 1.e-20 </gauss_seidel_tolerance>
</cmfd>

34
tests/regression_tests/cmfd_nofeed/results_true.dat
View file

@ -391,10 +391,6 @@ cmfd indices
1.000000E+00
1.000000E+00
k cmfd
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.170416E+00
1.172572E+00
1.171159E+00
@ -412,10 +408,6 @@ k cmfd
1.153714E+00
1.159485E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.203643E+00
3.204555E+00
3.210935E+00
@ -433,10 +425,6 @@ cmfd entropy
3.221580E+00
3.220523E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.009062E-03
4.869661E-03
2.997290E-03
@ -446,18 +434,14 @@ cmfd balance
1.403979E-03
1.398429E-03
1.818398E-03
1.761250E-03
1.646649E-03
1.480119E-03
1.399559E-03
1.400162E-03
1.178363E-03
1.292280E-03
1.458254E-03
1.437876E-03
1.276844E-03
1.289036E-03
1.330300E-03
1.139999E-03
1.233647E-03
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
5.397E-01
5.405E-01
5.412E-01
@ -475,10 +459,6 @@ cmfd dominance ratio
5.461E-01
5.443E-01
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.959834E-03
5.494667E-03
4.076255E-03

3
tests/regression_tests/cmfd_nofeed/settings.xml
View file

@ -23,7 +23,4 @@
</mesh>
<entropy_mesh>10</entropy_mesh>
<!-- Run CMFD -->
<run_cmfd>true</run_cmfd>
</settings>

22
tests/regression_tests/cmfd_nofeed/test.py
View file

@ -1,6 +1,26 @@
from tests.testing_harness import CMFDTestHarness
from openmc import cmfd
import numpy as np
def test_cmfd_nofeed():
harness = CMFDTestHarness('statepoint.20.h5')
"""Test 1 group CMFD solver without CMFD feedback"""
# Initialize and set CMFD mesh
cmfd_mesh = cmfd.CMFDMesh()
cmfd_mesh.lower_left = -10.0, -1.0, -1.0
cmfd_mesh.upper_right = 10.0, 1.0, 1.0
cmfd_mesh.dimension = 10, 1, 1
cmfd_mesh.albedo = 0.0, 0.0, 1.0, 1.0, 1.0, 1.0
# Initialize and run CMFDRun object
cmfd_run = cmfd.CMFDRun()
cmfd_run.mesh = cmfd_mesh
cmfd_run.begin = 5
cmfd_run.display = {'dominance': True}
cmfd_run.feedback = False
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
cmfd_run.run()
# Initialize and run CMFD test harness
harness = CMFDTestHarness('statepoint.20.h5', cmfd_run)
harness.main()

84
tests/testing_harness.py
View file

@ -133,35 +133,67 @@ class HashedTestHarness(TestHarness):
class CMFDTestHarness(TestHarness):
"""Specialized TestHarness for running OpenMC CMFD tests."""
def _get_results(self):
"""Digest info in the statepoint and return as a string."""
def __init__(self, statepoint_name, cmfd_run):
self._sp_name = statepoint_name
self._create_cmfd_result_str(cmfd_run)
# Write out the eigenvalue and tallies.
outstr = super()._get_results()
def _create_cmfd_result_str(self, cmfd_run):
"""Create CMFD result string from variables of CMFDRun instance"""
outstr = 'cmfd indices\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfd_run.indices])
outstr += '\nk cmfd\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfd_run.k_cmfd])
outstr += '\ncmfd entropy\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfd_run.entropy])
outstr += '\ncmfd balance\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfd_run.balance])
outstr += '\ncmfd dominance ratio\n'
outstr += '\n'.join(['{0:10.3E}'.format(x) for x in cmfd_run.dom])
outstr += '\ncmfd openmc source comparison\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfd_run.src_cmp])
outstr += '\ncmfd source\n'
cmfdsrc = np.reshape(cmfd_run.cmfd_src, np.product(cmfd_run.indices),
order='F')
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfdsrc])
outstr += '\n'
self._cmfdrun_results = outstr
# Read the statepoint file.
statepoint = glob.glob(self._sp_name)[0]
with openmc.StatePoint(statepoint) as sp:
# Write out CMFD data.
outstr += 'cmfd indices\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in sp.cmfd_indices])
outstr += '\nk cmfd\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in sp.k_cmfd])
outstr += '\ncmfd entropy\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in sp.cmfd_entropy])
outstr += '\ncmfd balance\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in sp.cmfd_balance])
outstr += '\ncmfd dominance ratio\n'
outstr += '\n'.join(['{0:10.3E}'.format(x) for x in sp.cmfd_dominance])
outstr += '\ncmfd openmc source comparison\n'
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in sp.cmfd_srccmp])
outstr += '\ncmfd source\n'
cmfdsrc = np.reshape(sp.cmfd_src, np.product(sp.cmfd_indices),
order='F')
outstr += '\n'.join(['{0:12.6E}'.format(x) for x in cmfdsrc])
outstr += '\n'
def execute_test(self):
"""Don't call _run_openmc as OpenMC will be called through C API for
CMFD tests, and write CMFD results that were passsed as argument
return outstr
"""
try:
self._test_output_created()
results = self._get_results()
results += self._cmfdrun_results
self._write_results(results)
self._compare_results()
finally:
self._cleanup()
def update_results(self):
"""Don't call _run_openmc as OpenMC will be called through C API for
CMFD tests, and write CMFD results that were passsed as argument
"""
try:
self._test_output_created()
results = self._get_results()
results += self._cmfdrun_results
self._write_results(results)
self._overwrite_results()
finally:
self._cleanup()
def _cleanup(self):
"""Delete output files for numpy matrices and flux vectors."""
super()._cleanup()
output = ['loss.npz', 'loss.dat', 'prod.npz', 'prod.dat',
'fluxvec.npy', 'fluxvec.dat']
for f in output:
if os.path.exists(f):
os.remove(f)
class ParticleRestartTestHarness(TestHarness):

1
tests/unit_tests/test_settings.py
View file

@ -21,7 +21,6 @@ def test_export_to_xml(run_in_tmpdir):
s.confidence_intervals = True
s.cross_sections = '/path/to/cross_sections.xml'
s.ptables = True
s.run_cmfd = False
s.seed = 17
s.survival_biasing = True
s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy': 1.0e-5}