From ee9a6ab065faff63150fca6d9aa3078ae8e57b62 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Sat, 24 Nov 2018 00:13:37 -0500 Subject: [PATCH] So begins the CMFD Fortran purge :D (~3.6K lines) --- CMakeLists.txt | 8 - docs/source/io_formats/cmfd.rst | 221 ------ docs/source/io_formats/index.rst | 1 - docs/source/io_formats/settings.rst | 10 - docs/source/io_formats/statepoint.rst | 20 - docs/source/pythonapi/base.rst | 63 +- docs/source/usersguide/basics.rst | 4 - openmc/cmfd.py | 515 +++----------- openmc/model/model.py | 18 +- openmc/settings.py | 18 - schemas.xml | 1 - src/api.F90 | 4 - src/cmfd_data.F90 | 943 -------------------------- src/cmfd_execute.F90 | 408 ----------- src/cmfd_execute.cpp | 43 -- src/cmfd_header.F90 | 272 -------- src/cmfd_input.F90 | 495 -------------- src/cmfd_loss_operator.F90 | 435 ------------ src/cmfd_prod_operator.F90 | 199 ------ src/cmfd_solver.F90 | 809 ---------------------- src/constants.F90 | 12 - src/input_xml.F90 | 4 - src/output.F90 | 44 -- src/relaxng/cmfd.rnc | 53 -- src/relaxng/cmfd.rng | 215 ------ src/relaxng/settings.rng | 5 - src/settings.F90 | 3 - src/simulation.cpp | 15 - src/state_point.F90 | 48 +- 29 files changed, 148 insertions(+), 4738 deletions(-) delete mode 100644 docs/source/io_formats/cmfd.rst delete mode 100644 src/cmfd_data.F90 delete mode 100644 src/cmfd_execute.F90 delete mode 100644 src/cmfd_execute.cpp delete mode 100644 src/cmfd_header.F90 delete mode 100644 src/cmfd_input.F90 delete mode 100644 src/cmfd_loss_operator.F90 delete mode 100644 src/cmfd_prod_operator.F90 delete mode 100644 src/cmfd_solver.F90 delete mode 100644 src/relaxng/cmfd.rnc delete mode 100644 src/relaxng/cmfd.rng diff --git a/CMakeLists.txt b/CMakeLists.txt index 7fc6fcec6..c9bb9f3c8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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,6 @@ add_library(libopenmc SHARED src/tallies/trigger_header.F90 src/dagmc.cpp src/cell.cpp - src/cmfd_execute.cpp src/cmfd_solver.cpp src/cross_sections.cpp src/distribution.cpp diff --git a/docs/source/io_formats/cmfd.rst b/docs/source/io_formats/cmfd.rst deleted file mode 100644 index 80e8365c7..000000000 --- a/docs/source/io_formats/cmfd.rst +++ /dev/null @@ -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 ```` element in -``settings.xml`` should be set to "true". - -------------------- -```` Element -------------------- - -The ```` element controls what batch CMFD calculations should begin. - - *Default*: 1 - ------------------------- -```` Element ------------------------- - -The ```` 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 - ---------------------- -```` Element ---------------------- - -The ```` 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 - -------------------------- -```` Element -------------------------- - -The ```` 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 - ----------------------- -```` Element ----------------------- - -The ```` 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 - ------------------------------------- -```` Element ------------------------------------- - -The ```` 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 - --------------------- -```` Element --------------------- - -The ```` element specifies the tolerance on the eigenvalue when performing -CMFD power iteration. - - *Default*: 1.e-8 - ------------------- -```` 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. - ------------------- -```` Element ------------------- - -The ```` 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 - ---------------------------- -```` Element ---------------------------- - -The ```` 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 - -------------------------- -```` Element -------------------------- - -The ```` element can be turned on with "true" to have an adjoint -calculation be performed on the last batch when CMFD is active. - - *Default*: false - --------------------- -```` Element --------------------- - -The ```` 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 - ----------------------- -```` Element ----------------------- - -The ```` 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 - ------------------- -```` Element ------------------- - -The ```` element specifies the tolerance on the fission source when performing -CMFD power iteration. - - *Default*: 1.e-8 - -------------------------- -```` Element -------------------------- - -The ```` element contains a list of batch numbers in which CMFD tallies -should be reset. - - *Default*: None - ----------------------------- -```` Element ----------------------------- - -The ```` 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 diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index ef83bd23d..7ae22c0a3 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -19,7 +19,6 @@ Input Files settings tallies plots - cmfd ---------- Data Files diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index f19a9e902..178c44080 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -361,16 +361,6 @@ or sub-elements: .. note:: This element is not used in the multi-group :ref:`energy_mode`. ----------------------- -```` Element ----------------------- - -The ```` 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 - ---------------------- ```` Element ---------------------- diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 38dae7dd7..9c5b8f526 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -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. diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 84bb3cc54..5202ea591 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -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,53 @@ 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 +------------------------------------------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.CMFDMesh + openmc.CMFDRun + +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. The following properties can be set through the CMFDRun class: + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myfunction.rst + + openmc.CMFDRun.cmfd_begin + openmc.CMFDRun.dhat_reset + openmc.CMFDRun.cmfd_display + openmc.CMFDRun.cmfd_downscatter + openmc.CMFDRun.cmfd_feedback + openmc.CMFDRun.cmfd_ktol + openmc.CMFDRun.cmfd_mesh + openmc.CMFDRun.norm + openmc.CMFDRun.cmfd_adjoint_type + openmc.CMFDRun.cmfd_power_monitor + openmc.CMFDRun.cmfd_run_adjoint + openmc.CMFDRun.cmfd_shift + openmc.CMFDRun.cmfd_stol + openmc.CMFDRun.cmfd_spectral + openmc.CMFDRun.cmfd_reset + openmc.CMFDRun.cmfd_write_matrices + openmc.CMFDRun.gauss_seidel_tolerance + +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 + diff --git a/docs/source/usersguide/basics.rst b/docs/source/usersguide/basics.rst index 74e65f771..2dd91ee0f 100644 --- a/docs/source/usersguide/basics.rst +++ b/docs/source/usersguide/basics.rst @@ -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) -------------------------------- diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 2c0a3eca2..ee9a8a1b8 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -12,7 +12,6 @@ References from collections.abc import Iterable from numbers import Real, Integral -from xml.etree import ElementTree as ET # TODO Remove import sys import numpy as np # Line below is added to suppress warnings when using numpy.divide to @@ -30,7 +29,6 @@ except ImportError: have_mpi = False import openmc.capi -from openmc._xml import clean_indentation #TODO Remove from openmc.checkvalue import (check_type, check_length, check_value, check_greater_than, check_less_than) from openmc.exceptions import OpenMCError @@ -242,385 +240,6 @@ class CMFDMesh(object): return element -# REMOVE this entire class -class CMFD(object): - r"""Parameters that control the use of coarse-mesh finite difference acceleration - in OpenMC. This corresponds directly to the cmfd.xml input file. - - Attributes - ---------- - begin : int - Batch number at which CMFD calculations should begin - dhat_reset : bool - Indicate whether :math:`\widehat{D}` nonlinear CMFD parameters should be - reset to zero before solving CMFD eigenproblem. - display : {'balance', 'dominance', 'entropy', 'source'} - Set 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. - * "entropy" - prints the *entropy* of the CMFD predicted fission source. - * "source" - prints the RMS [%] between the OpenMC fission source and - CMFD fission source. - downscatter : bool - Indicate whether an effective downscatter cross section should be used - when using 2-group CMFD. - feedback : bool - Indicate or not the CMFD diffusion result is used to adjust the weight - of fission source neutrons on the next OpenMC batch. Defaults to False. - gauss_seidel_tolerance : Iterable of float - Two parameters specifying the absolute inner tolerance and the relative - inner tolerance for Gauss-Seidel iterations when performing CMFD. - ktol : float - Tolerance on the eigenvalue when performing CMFD power iteration - cmfd_mesh : openmc.CMFDMesh - Structured mesh to be used for acceleration - norm : float - Normalization factor applied to the CMFD fission source distribution - power_monitor : bool - View convergence of power iteration during CMFD acceleration - run_adjoint : bool - Perform adjoint calculation on the last batch - shift : float - Optional Wielandt shift parameter for accelerating power iterations. By - default, it is very large so there is effectively no impact. - spectral : float - Optional spectral radius that can be used to accelerate the convergence - of Gauss-Seidel iterations during CMFD power iteration. - stol : float - Tolerance on the fission source when performing CMFD power iteration - tally_reset : list of int - List of batch numbers at which CMFD tallies should be reset - write_matrices : bool - Write sparse matrices that are used during CMFD acceleration (loss, - production) to file - - """ - - def __init__(self): - self._begin = None - self._dhat_reset = None - self._display = None - self._downscatter = None - self._feedback = None - self._gauss_seidel_tolerance = None - self._ktol = None - self._cmfd_mesh = None - self._norm = None - self._power_monitor = None - self._run_adjoint = None - self._shift = None - self._spectral = None - self._stol = None - self._tally_reset = None - self._write_matrices = None - - self._cmfd_file = ET.Element("cmfd") - self._cmfd_mesh_element = None - - @property - def begin(self): - return self._begin - - @property - def dhat_reset(self): - return self._dhat_reset - - @property - def display(self): - return self._display - - @property - def downscatter(self): - return self._downscatter - - @property - def feedback(self): - return self._feedback - - @property - def gauss_seidel_tolerance(self): - return self._gauss_seidel_tolerance - - @property - def ktol(self): - return self._ktol - - @property - def cmfd_mesh(self): - return self._cmfd_mesh - - @property - def norm(self): - return self._norm - - @property - def power_monitor(self): - return self._power_monitor - - @property - def run_adjoint(self): - return self._run_adjoint - - @property - def shift(self): - return self._shift - - @property - def spectral(self): - return self._spectral - - @property - def stol(self): - return self._stol - - @property - def tally_reset(self): - return self._tally_reset - - @property - def write_matrices(self): - return self._write_matrices - - @begin.setter - def begin(self, begin): - check_type('CMFD begin batch', begin, Integral) - check_greater_than('CMFD begin batch', begin, 0) - self._begin = begin - - @dhat_reset.setter - def dhat_reset(self, dhat_reset): - check_type('CMFD Dhat reset', dhat_reset, bool) - self._dhat_reset = dhat_reset - - @display.setter - def display(self, display): - check_type('CMFD display', display, str) - check_value('CMFD display', display, - ['balance', 'dominance', 'entropy', 'source']) - self._display = display - - @downscatter.setter - def downscatter(self, downscatter): - check_type('CMFD downscatter', downscatter, bool) - self._downscatter = downscatter - - @feedback.setter - def feedback(self, feedback): - check_type('CMFD feedback', feedback, bool) - self._feedback = feedback - - @gauss_seidel_tolerance.setter - def gauss_seidel_tolerance(self, gauss_seidel_tolerance): - check_type('CMFD Gauss-Seidel tolerance', gauss_seidel_tolerance, - Iterable, Real) - check_length('Gauss-Seidel tolerance', gauss_seidel_tolerance, 2) - self._gauss_seidel_tolerance = gauss_seidel_tolerance - - @ktol.setter - def ktol(self, ktol): - check_type('CMFD eigenvalue tolerance', ktol, Real) - self._ktol = ktol - - @cmfd_mesh.setter - def cmfd_mesh(self, mesh): - check_type('CMFD mesh', mesh, CMFDMesh) - - # Check dimension defined - if mesh.dimension is None: - raise ValueError('CMFD mesh requires spatial ' - 'dimensions to be specified') - - # Check lower left defined - if mesh.lower_left is None: - raise ValueError('CMFD mesh requires lower left coordinates ' - 'to be specified') - - # Check that both upper right and width both not defined - if mesh.upper_right is not None and mesh.width is not None: - raise ValueError('Both upper right coordinates and width ' - 'cannot be specified for CMFD mesh') - - # Check that at least one of width or upper right is defined - if mesh.upper_right is None and mesh.width is None: - raise ValueError('CMFD mesh requires either upper right ' - 'coordinates or width to be specified') - - # Check width and lower length are same dimension and define upper_right - if mesh.width is not None: - check_length('CMFD mesh width', mesh.width, len(mesh.lower_left)) - mesh.upper_right = np.array(mesh.lower_left) + \ - np.array(mesh.width) * np.array(mesh.dimension) - - # Check upper_right and lower length are same dimension and define width - elif mesh.upper_right is not None: - check_length('CMFD mesh upper right', mesh.upper_right, \ - len(mesh.lower_left)) - # Check upper right coordinates are greater than lower left - if np.any(np.array(mesh.upper_right) <= np.array(mesh.lower_left)): - raise ValueError('CMFD mesh requires upper right ' - 'coordinates to be greater than lower ' - 'left coordinates') - mesh.width = np.true_divide( - (np.array(mesh.upper_right) - np.array(mesh.lower_left)), \ - np.array(mesh.dimension)) - self._cmfd_mesh = mesh - - @norm.setter - def norm(self, norm): - check_type('CMFD norm', norm, Real) - self._norm = norm - - @power_monitor.setter - def power_monitor(self, power_monitor): - check_type('CMFD power monitor', power_monitor, bool) - self._power_monitor = power_monitor - - @run_adjoint.setter - def run_adjoint(self, run_adjoint): - check_type('CMFD run adjoint', run_adjoint, bool) - self._run_adjoint = run_adjoint - - @shift.setter - def shift(self, shift): - check_type('CMFD Wielandt shift', shift, Real) - self._shift = shift - - @spectral.setter - def spectral(self, spectral): - check_type('CMFD spectral radius', spectral, Real) - self._spectral = spectral - - @stol.setter - def stol(self, stol): - check_type('CMFD fission source tolerance', stol, Real) - self._stol = stol - - @tally_reset.setter - def tally_reset(self, tally_reset): - check_type('tally reset batches', tally_reset, Iterable, Integral) - self._tally_reset = tally_reset - - @write_matrices.setter - def write_matrices(self, write_matrices): - check_type('CMFD write matrices', write_matrices, bool) - self._write_matrices = write_matrices - - - def _create_begin_subelement(self): - if self._begin is not None: - element = ET.SubElement(self._cmfd_file, "begin") - element.text = str(self._begin) - - def _create_dhat_reset_subelement(self): - if self._dhat_reset is not None: - element = ET.SubElement(self._cmfd_file, "dhat_reset") - element.text = str(self._dhat_reset).lower() - - def _create_display_subelement(self): - if self._display is not None: - element = ET.SubElement(self._cmfd_file, "display") - element.text = str(self._display) - - def _create_downscatter_subelement(self): - if self._downscatter is not None: - element = ET.SubElement(self._cmfd_file, "downscatter") - element.text = str(self._downscatter).lower() - - def _create_feedback_subelement(self): - if self._feedback is not None: - element = ET.SubElement(self._cmfd_file, "feeback") - element.text = str(self._feedback).lower() - - def _create_gauss_seidel_tolerance_subelement(self): - if self._gauss_seidel_tolerance is not None: - element = ET.SubElement(self._cmfd_file, "gauss_seidel_tolerance") - element.text = ' '.join(map(str, self._gauss_seidel_tolerance)) - - def _create_ktol_subelement(self): - if self._ktol is not None: - element = ET.SubElement(self._ktol, "ktol") - element.text = str(self._ktol) - - def _create_mesh_subelement(self): - if self._cmfd_mesh is not None: - xml_element = self._cmfd_mesh._get_xml_element() - self._cmfd_file.append(xml_element) - - def _create_norm_subelement(self): - if self._norm is not None: - element = ET.SubElement(self._cmfd_file, "norm") - element.text = str(self._norm) - - def _create_power_monitor_subelement(self): - if self._power_monitor is not None: - element = ET.SubElement(self._cmfd_file, "power_monitor") - element.text = str(self._power_monitor).lower() - - def _create_run_adjoint_subelement(self): - if self._run_adjoint is not None: - element = ET.SubElement(self._cmfd_file, "run_adjoint") - element.text = str(self._run_adjoint).lower() - - def _create_shift_subelement(self): - if self._shift is not None: - element = ET.SubElement(self._shift, "shift") - element.text = str(self._shift) - - def _create_spectral_subelement(self): - if self._spectral is not None: - element = ET.SubElement(self._spectral, "spectral") - element.text = str(self._spectral) - - def _create_stol_subelement(self): - if self._stol is not None: - element = ET.SubElement(self._stol, "stol") - element.text = str(self._stol) - - def _create_tally_reset_subelement(self): - if self._tally_reset is not None: - element = ET.SubElement(self._tally_reset, "tally_reset") - element.text = ' '.join(map(str, self._tally_reset)) - - def _create_write_matrices_subelement(self): - if self._write_matrices is not None: - element = ET.SubElement(self._cmfd_file, "write_matrices") - element.text = str(self._write_matrices).lower() - - def export_to_xml(self): - """Create a cmfd.xml file using the class data that can be used for an OpenMC - simulation. - - """ - - self._create_begin_subelement() - self._create_dhat_reset_subelement() - self._create_display_subelement() - self._create_downscatter_subelement() - self._create_feedback_subelement() - self._create_gauss_seidel_tolerance_subelement() - self._create_ktol_subelement() - self._create_mesh_subelement() - self._create_norm_subelement() - self._create_power_monitor_subelement() - self._create_run_adjoint_subelement() - self._create_shift_subelement() - self._create_spectral_subelement() - self._create_stol_subelement() - self._create_tally_reset_subelement() - self._create_write_matrices_subelement() - - # Clean the indentation in the file to be user-readable - clean_indentation(self._cmfd_file) - - # Write the XML Tree to the cmfd.xml file - tree = ET.ElementTree(self._cmfd_file) - tree.write("cmfd.xml", xml_declaration=True, - encoding='utf-8', method="xml") - - class CMFDRun(object): r"""Class to run openmc with CMFD acceleration through the C API. Running openmc in this manner obviates the need for defining CMFD parameters @@ -629,27 +248,6 @@ class CMFDRun(object): Attributes ---------- - cmfd_begin : int - Batch number at which CMFD calculations should begin - dhat_reset : bool - Indicate whether :math:`\widehat{D}` nonlinear CMFD parameters should be - reset to zero before solving CMFD eigenproblem. - cmfd_display : {'balance', 'dominance', 'entropy', 'source'} - Set 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. - * "entropy" - prints the *entropy* of the CMFD predicted fission source. - * "source" - prints the RMS [%] between the OpenMC fission source and - CMFD fission source. - cmfd_downscatter : bool - Indicate whether an effective downscatter cross section should be used - when using 2-group CMFD. - cmfd_feedback : bool - Indicate or not the CMFD diffusion result is used to adjust the weight - of fission source neutrons on the next OpenMC batch. Defaults to False. cmfd_ktol : float Tolerance on the eigenvalue when performing CMFD power iteration cmfd_mesh : openmc.CMFDMesh @@ -668,9 +266,14 @@ class CMFDRun(object): cmfd_reset : list of int List of batch numbers at which CMFD tallies should be reset cmfd_write_matrices : bool - # TODO update this to read "resultant normalized flux vector" Write sparse matrices that are used during CMFD acceleration (loss, - production) and resultant flux vector phi to file + production) and resultant normalized flux vector phi to file + cmfd_spectral : float + Optional spectral radius that can be used to accelerate the convergence + of Gauss-Seidel iterations during CMFD power iteration. + gauss_seidel_tolerance : Iterable of float + Two parameters specifying the absolute inner tolerance and the relative + inner tolerance for Gauss-Seidel iterations when performing CMFD. indices : numpy.ndarray Stores spatial and group dimensions as [nx, ny, nz, ng] egrid : numpy.ndarray @@ -774,9 +377,6 @@ class CMFDRun(object): k_cmfd : list of floats List of CMFD k-effectives, stored for each generation that CMFD is invoked - spectral : float - Optional spectral radius that can be used to accelerate the convergence - of Gauss-Seidel iterations during CMFD power iteration. resnb : numpy.ndarray Residual from solving neutron balance equations time_cmfd : float @@ -787,13 +387,89 @@ class CMFDRun(object): Time for solving CMFD matrix equations, in seconds intracomm : mpi4py.MPI.Intracomm or None MPI intercommunicator for running MPI commands - - TODO Remove CMFD constants, timing variables in Fortran - TODO Remove cmfd_data.F90, cmfd_execute.F90, cmfd_header.F90, - cmfd_input.F90, cmfd_loss_operator.F90, cmfd_prod_operator.F90, - cmfd_solver.F90 - TODO Remove instances of cmfd in sourcepoint.F90, output.F90, - simulation.F90, api.F90, settings.F90, input_xml.F90, plots.F90 + first_x_accel : tuple + Indices in CMFD problem where first x element is an accelerated region + Precomputed and stored for updating CMFD arrays + last_x_accel : tuple + Indices in CMFD problem where last x element is an accelerated region + Precomputed and stored for updating CMFD arrays + first_y_accel : tuple + Indices in CMFD problem where first y element is an accelerated region + Precomputed and stored for updating CMFD arrays + last_y_accel : tuple + Indices in CMFD problem where last y element is an accelerated region + Precomputed and stored for updating CMFD arrays + first_z_accel : tuple + Indices in CMFD problem where first z element is an accelerated region + Precomputed and stored for updating CMFD arrays + last_z_accel : tuple + Indices in CMFD problem where last z element is an accelerated region + Precomputed and stored for updating CMFD arrays + notfirst_x_accel : tuple + Indices in CMFD problem where all x element excluding first are + accelerated regions. Precomputed and stored for updating CMFD arrays + notlast_x_accel : tuple + Indices in CMFD problem where all x element excluding last are + accelerated regions. Precomputed and stored for updating CMFD arrays + notfirst_y_accel : tuple + Indices in CMFD problem where all y element excluding first are + accelerated regions. Precomputed and stored for updating CMFD arrays + notlast_y_accel : tuple + Indices in CMFD problem where all y element excluding last are + accelerated regions. Precomputed and stored for updating CMFD arrays + notfirst_z_accel : tuple + Indices in CMFD problem where all z element excluding first are + accelerated regions. Precomputed and stored for updating CMFD arrays + notlast_z_accel : tuple + Indices in CMFD problem where all z element excluding last are + accelerated regions. Precomputed and stored for updating CMFD arrays + is_adj_ref_left : numpy.ndarray + Boolean array of all indices in notfirst_x_accel that neighbor a reflector + region to the left. Precomputed and stored for updating CMFD arrays + is_adj_ref_right : numpy.ndarray + Boolean array of all indices in notlast_x_accel that neighbor a reflector + region to the right. Precomputed and stored for updating CMFD arrays + is_adj_ref_back : numpy.ndarray + Boolean array of all indices in notfirst_y_accel that neighbor a reflector + region to the back. Precomputed and stored for updating CMFD arrays + is_adj_ref_front : numpy.ndarray + Boolean array of all indices in notlast_y_accel that neighbor a reflector + region to the front. Precomputed and stored for updating CMFD arrays + is_adj_ref_bottom : numpy.ndarray + Boolean array of all indices in notfirst_z_accel that neighbor a reflector + region to the bottom. Precomputed and stored for updating CMFD arrays + is_adj_ref_top : numpy.ndarray + Boolean array of all indices in notlast_z_accel that neighbor a reflector + region to the top. Precomputed and stored for updating CMFD arrays + accel_idxs : tuple + All indices in CMFD problem that are accelerated. Precomputed and + stored for updating CMFD matrixes + accel_neig_left_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the left + accel_neig_right_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the right + accel_neig_back_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the back + accel_neig_front_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the front + accel_neig_bottom_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the bottom + accel_neig_top_idxs : tuple + All indices in CMFD problem that are accelerated and have a neighbor to + the top + loss_row : numpy.ndarray + All row indices in loss matrix that have nonzero elements + loss_col : numpy.ndarray + All column indices in loss matrix that have nonzero elements + prod_row : numpy.ndarray + All row indices in production matrix that have nonzero elements + prod_col : numpy.ndarray + All column indices in production matrix that have nonzero elements """ @@ -860,8 +536,10 @@ class CMFDRun(object): self._time_cmfd = None self._time_cmfdbuild = None self._time_cmfdsolve = None + self._intracomm = None # Add all index-related variables, for numpy vectorization + self._first_x_accel = None self._last_x_accel = None self._first_y_accel = None self._last_y_accel = None @@ -890,7 +568,6 @@ class CMFDRun(object): self._loss_col = None self._prod_row = None self._prod_col = None - self._intracomm = None @property @@ -1440,11 +1117,9 @@ class CMFDRun(object): # Write out flux vector if self._cmfd_write_matrices: if adjoint: - # TODO Change to self._adj_phi - self._write_vector(phi, 'adj_fluxvec') + self._write_vector(self._adj_phi, 'adj_fluxvec') else: - # TODO Change to self._phi - self._write_vector(phi, 'fluxvec') + self._write_vector(self._phi, 'fluxvec') def _write_vector(self, vector, base_filename): """Write a 1-D numpy array to file and also save it in .npy format. This @@ -2364,7 +2039,7 @@ class CMFDRun(object): is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT x_inds, y_inds, z_inds = np.indices((nx, ny, nz)) - # Define slice equivalent to _accel[0,:,:] + # Define slice equivalent to is_accel[0,:,:] slice_x = x_inds[:1,:,:] slice_y = y_inds[:1,:,:] slice_z = z_inds[:1,:,:] diff --git a/openmc/model/model.py b/openmc/model/model.py index 8583da27a..28d19713c 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -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() diff --git a/openmc/settings.py b/openmc/settings.py index 2bbe0864e..1d97250ed 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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) diff --git a/schemas.xml b/schemas.xml index 465b9b398..3e586ec6a 100644 --- a/schemas.xml +++ b/schemas.xml @@ -5,6 +5,5 @@ - diff --git a/src/api.F90 b/src/api.F90 index efcfa16a6..30310e5fc 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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 @@ -103,9 +102,6 @@ contains call free_memory_dagmc() #endif - ! Deallocate CMFD - call deallocate_cmfd(cmfd) - end subroutine free_memory end module openmc_api diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 deleted file mode 100644 index f2fb4dc26..000000000 --- a/src/cmfd_data.F90 +++ /dev/null @@ -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 diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 deleted file mode 100644 index 681d08c5c..000000000 --- a/src/cmfd_execute.F90 +++ /dev/null @@ -1,408 +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, only: source_bank - 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 :: 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 - 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, int(work,4) - - ! Determine spatial bin - call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh) - - ! Determine energy bin - n_groups = size(cmfd % egrid) - 1 - if (source_bank(i) % E < cmfd % egrid(1)) then - e_bin = 1 - if (master) call warning('Source pt below energy grid') - elseif (source_bank(i) % E > 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(i) % E) - 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 - source_bank(i) % wgt = source_bank(i) % wgt * & - 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 diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp deleted file mode 100644 index 0af584c1b..000000000 --- a/src/cmfd_execute.cpp +++ /dev/null @@ -1,43 +0,0 @@ -#include // for copy -#include -#include - -#include "xtensor/xarray.hpp" -#include "xtensor/xio.hpp" - -#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 pointer to source bank - Bank* source_bank; - int64_t n; - openmc_source_bank(&source_bank, &n); - - // Get source counts in each mesh bin / energy bin - auto& m = model::meshes.at(settings::index_cmfd_mesh); - xt::xarray counts = m->count_sites(simulation::work, source_bank, 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 diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 deleted file mode 100644 index f4b14b6f2..000000000 --- a/src/cmfd_header.F90 +++ /dev/null @@ -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 diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 deleted file mode 100644 index d1e0f6a9a..000000000 --- a/src/cmfd_input.F90 +++ /dev/null @@ -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 diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 deleted file mode 100644 index 50a75e52d..000000000 --- a/src/cmfd_loss_operator.F90 +++ /dev/null @@ -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 diff --git a/src/cmfd_prod_operator.F90 b/src/cmfd_prod_operator.F90 deleted file mode 100644 index 91ce43b35..000000000 --- a/src/cmfd_prod_operator.F90 +++ /dev/null @@ -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 diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 deleted file mode 100644 index 01404d9b0..000000000 --- a/src/cmfd_solver.F90 +++ /dev/null @@ -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 diff --git a/src/constants.F90 b/src/constants.F90 index a9812e468..0e667c32a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1da1bf9df..32fd9f2a1 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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() diff --git a/src/output.F90 b/src/output.F90 index 485d957a6..d2130a87a 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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() diff --git a/src/relaxng/cmfd.rnc b/src/relaxng/cmfd.rnc deleted file mode 100644 index dde97f971..000000000 --- a/src/relaxng/cmfd.rnc +++ /dev/null @@ -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+ } }? - -} diff --git a/src/relaxng/cmfd.rng b/src/relaxng/cmfd.rng deleted file mode 100644 index dc096ab87..000000000 --- a/src/relaxng/cmfd.rng +++ /dev/null @@ -1,215 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 92a1d748a..3ed15c3ed 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -259,11 +259,6 @@ - - - - - diff --git a/src/settings.F90 b/src/settings.F90 index 6c8a61036..16f35341c 100644 --- a/src/settings.F90 +++ b/src/settings.F90 @@ -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 diff --git a/src/simulation.cpp b/src/simulation.cpp index 50bf12521..4a4f632eb 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -21,15 +21,10 @@ namespace openmc { // data/functions from Fortran side -extern "C" bool cmfd_on; - extern "C" void accumulate_tallies(); extern "C" void allocate_banks(); 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 join_bank_from_threads(); extern "C" void load_state_point(); @@ -94,9 +89,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(); @@ -316,11 +308,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(); } @@ -340,8 +327,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(); } diff --git a/src/state_point.F90 b/src/state_point.F90 index 2692856dc..0d33d52d1 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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 @@ -64,7 +63,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(:) @@ -169,23 +168,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") @@ -413,13 +395,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) @@ -453,7 +428,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 @@ -548,26 +522,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