diff --git a/examples/cmfd_testing/basic/results_true.dat b/examples/cmfd_testing/basic/results_true.dat deleted file mode 100644 index aba219ba8..000000000 --- a/examples/cmfd_testing/basic/results_true.dat +++ /dev/null @@ -1,508 +0,0 @@ -k-combined: -1.169891E+00 6.289489E-03 -tally 1: -1.173922E+01 -1.385461E+01 -2.164076E+01 -4.699368E+01 -2.906462E+01 -8.464937E+01 -3.382312E+01 -1.147095E+02 -3.632006E+01 -1.323878E+02 -3.655413E+01 -1.341064E+02 -3.347757E+01 -1.124264E+02 -2.931336E+01 -8.607239E+01 -2.182947E+01 -4.789565E+01 -1.147668E+01 -1.325716E+01 -tally 2: -2.298190E+01 -2.667071E+01 -1.600292E+01 -1.293670E+01 -4.268506E+01 -9.161216E+01 -3.022909E+01 -4.598915E+01 -5.680399E+01 -1.623879E+02 -4.033805E+01 -8.196263E+01 -6.814742E+01 -2.331778E+02 -4.851618E+01 -1.182330E+02 -7.392923E+01 -2.740255E+02 -5.253586E+01 -1.384152E+02 -7.332860E+01 -2.698608E+02 -5.227405E+01 -1.371810E+02 -6.830172E+01 -2.340687E+02 -4.867159E+01 -1.188724E+02 -5.885634E+01 -1.736180E+02 -4.170434E+01 -8.719622E+01 -4.371848E+01 -9.592893E+01 -3.106403E+01 -4.844308E+01 -2.338413E+01 -2.752467E+01 -1.636713E+01 -1.347770E+01 -tally 3: -1.538752E+01 -1.196478E+01 -1.079685E+00 -6.010786E-02 -2.911906E+01 -4.269070E+01 -1.822657E+00 -1.671850E-01 -3.885421E+01 -7.608218E+01 -2.541516E+00 -3.262451E-01 -4.673300E+01 -1.097036E+02 -2.885307E+00 -4.214444E-01 -5.059247E+01 -1.283984E+02 -3.222796E+00 -5.237329E-01 -5.034856E+01 -1.272538E+02 -3.230225E+00 -5.273424E-01 -4.688476E+01 -1.103152E+02 -2.941287E+00 -4.363749E-01 -4.013746E+01 -8.077506E+01 -2.634234E+00 -3.520270E-01 -2.996887E+01 -4.509953E+01 -1.946504E+00 -1.919104E-01 -1.575260E+01 -1.248707E+01 -1.020705E+00 -5.413569E-02 -tally 4: -3.049469E+00 -4.677325E-01 -0.000000E+00 -0.000000E+00 -2.770358E+00 -3.879191E-01 -5.514939E+00 -1.528899E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -5.514939E+00 -1.528899E+00 -2.770358E+00 -3.879191E-01 -5.032131E+00 -1.275040E+00 -7.294002E+00 -2.675589E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -7.294002E+00 -2.675589E+00 -5.032131E+00 -1.275040E+00 -7.036008E+00 -2.490718E+00 -8.668860E+00 -3.776102E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -8.668860E+00 -3.776102E+00 -7.036008E+00 -2.490718E+00 -8.352414E+00 -3.501945E+00 -9.345868E+00 -4.380719E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -9.345868E+00 -4.380719E+00 -8.352414E+00 -3.501945E+00 -9.093766E+00 -4.158282E+00 -9.223771E+00 -4.270120E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -9.223771E+00 -4.270120E+00 -9.093766E+00 -4.158282E+00 -9.219150E+00 -4.264346E+00 -8.530966E+00 -3.651778E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -8.530966E+00 -3.651778E+00 -9.219150E+00 -4.264346E+00 -8.690373E+00 -3.785262E+00 -7.204424E+00 -2.604203E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -7.204424E+00 -2.604203E+00 -8.690373E+00 -3.785262E+00 -7.513640E+00 -2.833028E+00 -5.326721E+00 -1.426975E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -5.326721E+00 -1.426975E+00 -7.513640E+00 -2.833028E+00 -5.662215E+00 -1.607757E+00 -2.848381E+00 -4.093396E-01 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -2.848381E+00 -4.093396E-01 -5.662215E+00 -1.607757E+00 -3.025812E+00 -4.597241E-01 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -tally 5: -1.538652E+01 -1.196332E+01 -2.252427E+00 -2.605738E-01 -2.911344E+01 -4.267319E+01 -3.873926E+00 -7.615035E-01 -3.884516E+01 -7.604619E+01 -5.280610E+00 -1.414008E+00 -4.672391E+01 -1.096625E+02 -6.261805E+00 -1.983205E+00 -5.058447E+01 -1.283588E+02 -6.733810E+00 -2.278242E+00 -5.033589E+01 -1.271898E+02 -6.714658E+00 -2.273652E+00 -4.687563E+01 -1.102719E+02 -6.215002E+00 -1.956978E+00 -4.013134E+01 -8.075062E+01 -5.253064E+00 -1.396224E+00 -2.996497E+01 -4.508840E+01 -3.818076E+00 -7.509442E-01 -1.574994E+01 -1.248291E+01 -2.219928E+00 -2.515492E-01 -cmfd indices -1.000000E+01 -1.000000E+00 -1.000000E+00 -1.000000E+00 -k cmfd -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -1.170416E+00 -1.172966E+00 -1.165537E+00 -1.170979E+00 -1.161922E+00 -1.157523E+00 -1.158873E+00 -1.162877E+00 -1.167101E+00 -1.168130E+00 -1.170570E+00 -1.168115E+00 -1.174081E+00 -1.169458E+00 -1.167848E+00 -1.165116E+00 -cmfd entropy -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -3.203643E+00 -3.207943E+00 -3.213367E+00 -3.214360E+00 -3.219634E+00 -3.222232E+00 -3.221744E+00 -3.224544E+00 -3.225990E+00 -3.227769E+00 -3.227417E+00 -3.230728E+00 -3.231662E+00 -3.233316E+00 -3.233193E+00 -3.232564E+00 -cmfd balance -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -4.009062E-03 -4.431773E-03 -3.152666E-03 -3.510383E-03 -2.052089E-03 -2.068651E-03 -1.502427E-03 -1.589825E-03 -1.566020E-03 -1.219160E-03 -1.017888E-03 -9.771622E-04 -1.010120E-03 -1.073382E-03 -1.172758E-03 -9.827332E-04 -cmfd dominance ratio - 0.000E+00 - 0.000E+00 - 0.000E+00 - 0.000E+00 - 5.397E-01 - 5.425E-01 - 5.481E-01 - 5.473E-01 - 5.503E-01 - 5.502E-01 - 5.483E-01 - 5.520E-01 - 5.505E-01 - 3.216E-01 - 5.373E-01 - 5.517E-01 - 5.508E-01 - 5.524E-01 - 5.524E-01 - 5.523E-01 -cmfd openmc source comparison -0.000000E+00 -0.000000E+00 -0.000000E+00 -0.000000E+00 -6.959834E-03 -5.655657E-03 -3.886185E-03 -4.035116E-03 -3.043277E-03 -5.455475E-03 -4.515311E-03 -2.439840E-03 -2.114032E-03 -2.673132E-03 -2.431749E-03 -4.330928E-03 -3.404647E-03 -3.680298E-03 -3.309620E-03 -3.705541E-03 -cmfd source -4.697085E-02 -7.920706E-02 -1.107968E-01 -1.250932E-01 -1.383930E-01 -1.380648E-01 -1.246874E-01 -1.113705E-01 -8.203754E-02 -4.337882E-02 diff --git a/examples/cmfd_testing/basic/capi/cmfd_ref.xml b/examples/cmfd_testing/cmfd-feed/capi/cmfd_ref.xml similarity index 100% rename from examples/cmfd_testing/basic/capi/cmfd_ref.xml rename to examples/cmfd_testing/cmfd-feed/capi/cmfd_ref.xml diff --git a/examples/cmfd_testing/basic/capi/geometry.xml b/examples/cmfd_testing/cmfd-feed/capi/geometry.xml similarity index 100% rename from examples/cmfd_testing/basic/capi/geometry.xml rename to examples/cmfd_testing/cmfd-feed/capi/geometry.xml diff --git a/examples/cmfd_testing/basic/capi/materials.xml b/examples/cmfd_testing/cmfd-feed/capi/materials.xml similarity index 100% rename from examples/cmfd_testing/basic/capi/materials.xml rename to examples/cmfd_testing/cmfd-feed/capi/materials.xml diff --git a/examples/cmfd_testing/basic/capi/run_openmc_cmfd.py b/examples/cmfd_testing/cmfd-feed/capi/run_openmc_cmfd.py similarity index 68% rename from examples/cmfd_testing/basic/capi/run_openmc_cmfd.py rename to examples/cmfd_testing/cmfd-feed/capi/run_openmc_cmfd.py index 5db5394ac..73d8e9a59 100644 --- a/examples/cmfd_testing/basic/capi/run_openmc_cmfd.py +++ b/examples/cmfd_testing/cmfd-feed/capi/run_openmc_cmfd.py @@ -7,6 +7,8 @@ cmfd_mesh.lower_left = [-10, -1, -1] cmfd_mesh.upper_right = [10, 1, 1] cmfd_mesh.dimension = [10, 1, 1] cmfd_mesh.albedo = [0., 0., 1., 1., 1., 1.] +cmfd_mesh.energy = [0., 10000., 10000000.] +cmfd_mesh.map = [0,1,1,1,1,1,1,1,1,0] # Initialize CMFDRun object cmfd_run = openmc.CMFDRun() @@ -14,7 +16,10 @@ cmfd_run = openmc.CMFDRun() # Set all runtime parameters (cmfd_mesh, tolerances, tally_resets, etc) # All error checking done under the hood when setter function called cmfd_run.cmfd_mesh = cmfd_mesh -cmfd_run.cmfd_reset = [5,10] +cmfd_run.cmfd_begin = 5 +cmfd_run.cmfd_display = 'dominance' +cmfd_run.cmfd_feedback = True +cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] # Run CMFD cmfd_run.run() diff --git a/examples/cmfd_testing/basic/capi/settings.xml b/examples/cmfd_testing/cmfd-feed/capi/settings.xml similarity index 100% rename from examples/cmfd_testing/basic/capi/settings.xml rename to examples/cmfd_testing/cmfd-feed/capi/settings.xml diff --git a/examples/cmfd_testing/basic/capi/tallies.xml b/examples/cmfd_testing/cmfd-feed/capi/tallies.xml similarity index 100% rename from examples/cmfd_testing/basic/capi/tallies.xml rename to examples/cmfd_testing/cmfd-feed/capi/tallies.xml diff --git a/examples/cmfd_testing/basic/trad/cmfd.xml b/examples/cmfd_testing/cmfd-feed/trad/cmfd.xml similarity index 75% rename from examples/cmfd_testing/basic/trad/cmfd.xml rename to examples/cmfd_testing/cmfd-feed/trad/cmfd.xml index b1c623ef2..adda40340 100644 --- a/examples/cmfd_testing/basic/trad/cmfd.xml +++ b/examples/cmfd_testing/cmfd-feed/trad/cmfd.xml @@ -6,10 +6,11 @@ 10 1 1 10 1 1 0.0 0.0 1.0 1.0 1.0 1.0 + 0.0 10000.0 10000000.0 + 1 2 2 2 2 2 2 2 2 1 - 5 10 5 dominance - true + false 1.e-15 1.e-20 diff --git a/examples/cmfd_testing/basic/trad/geometry.xml b/examples/cmfd_testing/cmfd-feed/trad/geometry.xml similarity index 100% rename from examples/cmfd_testing/basic/trad/geometry.xml rename to examples/cmfd_testing/cmfd-feed/trad/geometry.xml diff --git a/examples/cmfd_testing/basic/trad/materials.xml b/examples/cmfd_testing/cmfd-feed/trad/materials.xml similarity index 100% rename from examples/cmfd_testing/basic/trad/materials.xml rename to examples/cmfd_testing/cmfd-feed/trad/materials.xml diff --git a/examples/cmfd_testing/basic/trad/settings.xml b/examples/cmfd_testing/cmfd-feed/trad/settings.xml similarity index 100% rename from examples/cmfd_testing/basic/trad/settings.xml rename to examples/cmfd_testing/cmfd-feed/trad/settings.xml diff --git a/examples/cmfd_testing/basic/trad/tallies.xml b/examples/cmfd_testing/cmfd-feed/trad/tallies.xml similarity index 100% rename from examples/cmfd_testing/basic/trad/tallies.xml rename to examples/cmfd_testing/cmfd-feed/trad/tallies.xml diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index d706112c4..61f0ee221 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -1,4 +1,5 @@ -from ctypes import c_int, c_int32, c_int64, c_double, c_char_p, POINTER +from ctypes import (c_int, c_int32, c_int64, c_double, c_char_p, c_bool, + POINTER) from . import _dll from .core import _DLLGlobal @@ -17,9 +18,13 @@ _dll.openmc_get_seed.restype = c_int64 class _Settings(object): # Attributes that are accessed through a descriptor batches = _DLLGlobal(c_int32, 'n_batches') + current_batch = _DLLGlobal(c_int, 'openmc_current_batch') generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') + master = _DLLGlobal(c_bool, 'openmc_master') particles = _DLLGlobal(c_int64, 'n_particles') + restart_batch = _DLLGlobal(c_int, 'openmc_restart_batch') + restart_run = _DLLGlobal(c_bool, 'openmc_restart_run') verbosity = _DLLGlobal(c_int, 'openmc_verbosity') @property diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 3cf5e6cbe..51af38634 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -10,6 +10,7 @@ References """ +# TODO: Check to make sure no redundant import statements from collections.abc import Iterable from numbers import Real, Integral from xml.etree import ElementTree as ET @@ -21,9 +22,28 @@ from openmc.clean_xml import clean_xml_indentation from openmc.checkvalue import (check_type, check_length, check_value, check_greater_than, check_less_than) -# Maximum/minimum neutron energies, from src/api.F90 -ENERGY_MAX_NEUTRON = np.inf -ENERGY_MIN_NEUTRON = 0. + +""" +-------------- +CMFD CONSTANTS +-------------- +""" +# Maximum/minimum neutron energies +_ENERGY_MAX_NEUTRON = np.inf +_ENERGY_MIN_NEUTRON = 0. + +# Tolerance for detecting zero flux values +_TINY_BIT = 1.e-8 + +# For non-accelerated regions on coarse mesh overlay +_CMFD_NOACCEL = 99999 + +# Constant to represent a zero flux "albedo" +_ZERO_FLUX = 999.0 + +# Constant for writing out no residual +_CMFD_NORES = 99999.0 + class CMFDMesh(object): """A structured Cartesian mesh used for Coarse Mesh Finite Difference (CMFD) @@ -156,7 +176,7 @@ class CMFDMesh(object): def map(self, meshmap): check_type('CMFD mesh map', meshmap, Iterable, Integral) for m in meshmap: - check_value('CMFD mesh map', m, [1, 2]) + check_value('CMFD mesh map', m, [0, 1]) self._map = meshmap def _get_xml_element(self): @@ -186,7 +206,7 @@ class CMFDMesh(object): if self.map is not None: subelement = ET.SubElement(element, "map") - subelement.text = ' '.join(map(str, self.map)) + subelement.text = ' '.join(map(str, [self.map[i]+1 for i in range(len(self.map))])) return element @@ -409,8 +429,8 @@ class CMFD(object): # 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') + '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)) @@ -577,11 +597,9 @@ class CMFDRun(object): Attributes ---------- - To add: cmfd_coremap: Flag for active core map - indices: Stores spatial and group dimensions as [nx, ny, nz, ng] + To add: indices: Stores spatial and group dimensions as [nx, ny, nz, ng] egrid: energy grid used for CMFD acceleration albedo: Albedo for global boundary conditions, taken from CMFD mesh. Set to [1,1,1,1,1,1] if not specified by user - cmfd_coremap: Optional acceleration map to overlay on coarse mesh spatial grid, taken from CMFD mesh n_cmfd_resets: Number of elements in tally_reset, list that stores batches where CMFD tallies should be reset cmfd_atoli: Absolute GS tolerance, set by gauss_seidel_tolerance cmfd_rtoli: Relative GS tolerance, set by gauss_seidel_tolerance @@ -592,7 +610,6 @@ class CMFDRun(object): Set to true if user specifies energy grid in CMFDMesh, false otherwise cmfd_on: Boolean to tell if cmfd solver should be initiated, based on whether current batch has reached variable cmfd_begin - batch_num: Current batch # Look at cmfd_header.F90 for description flux totalxs @@ -613,9 +630,12 @@ class CMFDRun(object): src_cmp dom k_cmfd + keff_bal + mat_dim TODO: Put descriptions for all methods in CMFDRun TODO All timing variables + TODO Get rid of CMFD constants """ @@ -641,7 +661,6 @@ class CMFDRun(object): self._cmfd_write_matrices = False # External variables used during runtime but users don't have control over - self._cmfd_coremap = False self._indices = np.zeros(4, dtype=int) self._egrid = None self._albedo = None @@ -654,7 +673,8 @@ class CMFDRun(object): self._cmfd_tally_ids = None self._energy_filters = None self._cmfd_on = False - self._batch_num = 1 + self._mat_dim = _CMFD_NOACCEL + self._keff_bal = None # Numpy arrays used to build CMFD matrices self._flux = None @@ -745,7 +765,7 @@ class CMFDRun(object): def cmfd_begin(self, cmfd_begin): check_type('CMFD begin batch', cmfd_begin, Integral) check_greater_than('CMFD begin batch', cmfd_begin, 0) - self._cmfd_begin = begin + self._cmfd_begin = cmfd_begin @dhat_reset.setter def dhat_reset(self, dhat_reset): @@ -793,12 +813,12 @@ class CMFDRun(object): # Check lower left defined if mesh.lower_left is None: raise ValueError('CMFD mesh requires lower left coordinates ' - 'to be specified') + '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') + '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: @@ -865,7 +885,8 @@ class CMFDRun(object): check_type('CMFD write matrices', cmfd_write_matrices, bool) self._cmfd_write_matrices = cmfd_write_matrices - def run(self): + def run(self, mpi_procs=None, omp_threads=None): + # TODO: Add logic for mpi_procs, omp_threads as args openmc.capi.init() self._configure_cmfd() openmc.capi.simulation_init() @@ -877,13 +898,13 @@ class CMFDRun(object): # CMFD update or skip it entirely if it is a restart run status = openmc.capi.next_batch_between_cmfd_init_execute() if status != 0: - self._cmfd_execute() + if self._cmfd_on: + self._execute_cmfd() # Status now determines whether another batch should be run # or simulation should be terminated. status = openmc.capi.next_batch_after_cmfd_execute() if status != 0: break - self._batch_num += 1 openmc.capi.simulation_finalize() openmc.capi.finalize() @@ -902,12 +923,14 @@ class CMFDRun(object): self._allocate_cmfd() def _read_cmfd_input(self): + # TODO: Print message with verbosity + # Print message print(' Configuring CMFD parameters for simulation') # Check if CMFD mesh is defined if self._cmfd_mesh is None: raise ValueError('No CMFD mesh has been specified for ' - 'simulation') + 'simulation') # Set spatial dimensions of CMFD object # Iterate through each element of self._cmfd_mesh.dimension @@ -923,7 +946,7 @@ class CMFDRun(object): self._energy_filters = True # TODO: MG mode check else: - self._egrid = np.array([ENERGY_MIN_NEUTRON, ENERGY_MAX_NEUTRON]) + self._egrid = np.array([_ENERGY_MIN_NEUTRON, _ENERGY_MAX_NEUTRON]) self._indices[3] = 1 self._energy_filters = False @@ -933,14 +956,16 @@ class CMFDRun(object): else: self._albedo = np.array([1.,1.,1.,1.,1.,1.]) - # Get acceleration map + # Get acceleration map, otherwise set all regions to be accelerated if self._cmfd_mesh.map is not None: check_length('CMFD coremap', self._cmfd_mesh.map, np.product(self._indices[0:3])) self._coremap = np.array(self._cmfd_mesh.map).reshape(( \ self._indices[0], self._indices[1], \ self._indices[2])) - self._cmfd_coremap = True + else: + self._coremap = np.ones((self._indices[0], self._indices[1], \ + self._indices[2])) # Set number of batches where cmfd tallies should be reset if self._cmfd_reset is not None: @@ -954,9 +979,6 @@ class CMFDRun(object): self._create_cmfd_tally() def _allocate_cmfd(self): - # Determine number of batches through C API - n_batches = openmc.capi.settings.batches - # Extract spatial and energy indices nx = self._indices[0] ny = self._indices[1] @@ -964,57 +986,92 @@ class CMFDRun(object): ng = self._indices[3] # Allocate flux, cross sections and diffusion coefficient - self._flux = np.zeros((ng, nx, ny, nz)) - self._totalxs = np.zeros((ng, nx, ny, nz)) - self._p1scattxs = np.zeros((ng, nx, ny, nz)) - self._scattxs = np.zeros((ng, ng, nx, ny, nz)) - self._nfissxs = np.zeros((ng, ng, nx, ny, nz)) - self._diffcof = np.zeros((ng, nx, ny, nz)) + self._flux = np.zeros((nx, ny, nz, ng)) + self._totalxs = np.zeros((nx, ny, nz, ng)) + self._p1scattxs = np.zeros((nx, ny, nz, ng)) + self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Outgoing, incoming + self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Outgoing, incoming + self._diffcof = np.zeros((nx, ny, nz, ng)) # Allocate dtilde and dhat - self._dtilde = np.zeros((6, ng, nx, ny, nz)) - self._dhat = np.zeros((6, ng, nx, ny, nz)) + self._dtilde = np.zeros((6, nx, ny, nz, ng)) + self._dhat = np.zeros((6, nx, ny, nz, ng)) # Allocate dimensions for each box (here for general case) self._hxyz = np.zeros((3, nx, ny, nz)) # Allocate surface currents - self._current = np.zeros((12, ng, nx, ny, nz)) + self._current = np.zeros((nx, ny, nz, ng, 12)) # Allocate source distributions - self._cmfd_src = np.zeros((ng, nx, ny, nz)) - self._openmc_src = np.zeros((ng, nx, ny, nz)) + self._cmfd_src = np.zeros((nx, ny, nz, ng)) + self._openmc_src = np.zeros((nx, ny, nz, ng)) # Allocate source weight modification variables - self._sourcecounts = np.zeros((ng, nx*ny*nz)) - self._weightfactors = np.ones((ng, nx, ny, nz)) + self._sourcecounts = np.zeros((nx*ny*nz, ng)) + self._weightfactors = np.ones((nx, ny, nz, ng)) # Allocate batchwise parameters - self._entropy = np.zeros((n_batches, )) - self._balance = np.zeros((n_batches, )) - self._src_cmp = np.zeros((n_batches, )) - self._dom = np.zeros((n_batches, )) - self._k_cmfd = np.zeros((n_batches, )) + self._entropy = [] + self._balance = [] + self._src_cmp = [] + self._dom = [] + self._k_cmfd = [] def _cmfd_init_batch(self): - if self._cmfd_begin == self._batch_num: + current_batch = openmc.capi.settings.current_batch + restart_run = openmc.capi.settings.restart_run + restart_batch = openmc.capi.settings.restart_batch + + if self._cmfd_begin == current_batch: self._cmfd_on = True - # TODO: Figure out how to tell if part of restart run - # if restart_run: - # return - # If this is a restart run and we are just replaying batches leave - # if (restart_run .and. current_batch <= restart_batch) return + # TODO: Test restart_batch + if restart_run and current_batch <= restart_batch: + return # Check to reset tallies - if self._n_cmfd_resets > 0 and self._batch_num in self._cmfd_reset: + if (self._n_cmfd_resets > 0 + and current_batch in self._cmfd_reset): self._cmfd_tally_reset() - def _cmfd_execute(self): - pass + def _execute_cmfd(self): + # CMFD single processor on master + if openmc.capi.settings.master: + + # TODO + #! Start cmfd timer + #call time_cmfd % start() + + # Create cmfd data from OpenMC tallies + self._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() + ''' def _cmfd_tally_reset(self): - # TODO: How does write_message in error.F90 work? + # TODO: Print message with verbosity # Print message print(' CMFD tallies reset') @@ -1023,6 +1080,168 @@ class CMFDRun(object): for tally_id in self._cmfd_tally_ids: tallies[tally_id].reset() + def _set_up_cmfd(self): + # Check for core map and set it up + if (self._mat_dim == _CMFD_NOACCEL): + self._set_coremap() + + # Calculate all cross sections based on reaction rates from last batch + self._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() + ''' + + def _set_coremap(self): + self._mat_dim = np.sum(self._coremap) + self._coremap = np.where(self._coremap == 0, + _CMFD_NOACCEL, self._coremap) + + def _compute_xs(self): + # Extract energy indices + ng = self._indices[3] + + # Set flux object and source distribution all to zeros + self._flux.fill(0.) + self._openmc_src.fill(0.) + + # Reset keff_bal to zero + self._keff_bal = 0. + + # Get tallies in-memory + tallies = openmc.capi.tallies + + # Set conditional numpy array as boolean vector based on coremap + # Repeat each value for number of groups in problem + is_cmfd_accel = np.repeat(self._coremap.ravel() != _CMFD_NOACCEL, ng) + + # Get flux from CMFD tally 0 + tally_id = self._cmfd_tally_ids[0] + tally_results = tallies[tally_id].results[:,0,1] + flux = np.where(is_cmfd_accel, tally_results, 0.) + + # Detect zero flux, abort if located + if np.any(flux[is_cmfd_accel] < _TINY_BIT): + # Get index of zero flux in flux array + idx = np.argmax(np.where(is_cmfd_accel, flux, 1) < _TINY_BIT) + + # Convert scalar idx to index in flux matrix + mat_idx = np.unravel_index(idx, self._flux.shape) + + # Throw error message (one-based indexing) + # Index of group is flipped + err_message = 'Detected zero flux without coremap overlay' + \ + ' at mesh: (' + \ + ', '.join(str(i+1) for i in mat_idx[:-1]) + \ + ') in group ' + str(ng-mat_idx[-1]) + raise ValueError(err_message) + + # Store flux and reshape + # Flux is flipped in energy axis as tally results are given in reverse + # order of energy group + self._flux = np.flip(flux.reshape(self._flux.shape), axis=3) + + # Get total rr and convert to total xs from CMFD tally 0 + tally_results = tallies[tally_id].results[:,1,1] + totalxs = np.divide(tally_results, flux, \ + where=flux>0, out=np.zeros_like(tally_results)) + + # Store total xs and reshape + # Total xs is flipped in energy axis as tally results are given in + # reverse order of energy group + self._totalxs = np.flip(totalxs.reshape(self._totalxs.shape), axis=3) + + # Get scattering xs from CMFD tally 1 + # flux is repeated to account for extra dimensionality of scattering xs + tally_id = self._cmfd_tally_ids[1] + tally_results = tallies[tally_id].results[:,0,1] + scattxs = np.divide(tally_results, \ + np.repeat(flux, ng), \ + where=np.repeat(flux>0, ng), \ + out=np.zeros_like(tally_results)) + + # Store scattxs and reshape + # Scattering xs is flipped in both incoming and outgoing energy axes + # as tally results are given in reverse order of energy group + self._scattxs = np.flip(scattxs.reshape(self._scattxs.shape), axis=3) + self._scattxs = np.flip(self._scattxs.reshape(self._scattxs.shape), \ + axis=4) + + # Get nu-fission xs from CMFD tally 1 + # flux is repeated to account for extra dimensionality of nu-fission xs + tally_results = tallies[tally_id].results[:,1,1] + num_realizations = tallies[tally_id].num_realizations + nfissxs = np.divide(tally_results, \ + np.repeat(flux, ng), \ + where=np.repeat(flux>0, ng), \ + out=np.zeros_like(tally_results)) + + # Store nfissxs and reshape + # Nu-fission xs is flipped in both incoming and outgoing energy axes + # as tally results are given in reverse order of energy group + self._nfissxs = np.flip(nfissxs.reshape(self._nfissxs.shape), axis=3) + self._nfissxs = np.flip(self._nfissxs.reshape(self._nfissxs.shape), \ + axis=4) + + # Filter nu-fission tally results to compute openmc source distribution + tally_results = np.where(np.repeat(flux>0, ng), tally_results, \ + 0.) + + # Openmc source distribution is sum of nu-fission rr in outgoing energies + openmc_src = np.sum(tally_results.reshape(self._nfissxs.shape), + axis=3) + + # Store openmc_src + # Openmc source is flipped in energy axis as tally results are given + # in reverse order of energy group + self._openmc_src = np.flip(openmc_src, axis=3) + + # Compute k_eff from source distribution + self._keff_bal = np.sum(self._openmc_src) / num_realizations + + # Normalize openmc source distribution + self._openmc_src /= np.sum(self._openmc_src) * self._norm + + # Get surface currents from CMFD tally 2 + tally_id = self._cmfd_tally_ids[2] + tally_results = tallies[tally_id].results[:,0,1] + + # Filter tally results to include only accelerated regions + tally_results = np.where(np.repeat(flux>0, 12), tally_results, 0.) + + # Reshape and store current + # Current is flipped in energy axis as tally results are given in + # reverse order of energy group + self._current = np.flip(tally_results.reshape(self._current.shape), \ + axis=3) + + # Get p1 scatter xs from CMFD tally 3 + tally_id = self._cmfd_tally_ids[3] + tally_results = tallies[tally_id].results[:,0,1] + + # Reshape and extract only p1 data from tally results (no need for p0 data) + p1scattrr = tally_results.reshape(self._p1scattxs.shape+(2,))[:,:,:,1] + + # Store p1 scatter xs + # p1 scatter xs is flipped in energy axis as tally results are given in + # reverse order of energy group + self._p1scattxs = np.divide(np.flip(p1scattrr, axis=3), self._flux, \ + where=self._flux>0, \ + out=np.zeros_like(p1scattrr)) + + # Calculate and store diffusion coefficient + self._diffcof = np.where(self._flux > 0, 1.0 / (3.0 * \ + (self._totalxs - self._p1scattxs)), 0.) + def _create_cmfd_tally(self): # Create Mesh object based on CMFDMesh, stored internally cmfd_mesh = openmc.capi.Mesh() @@ -1083,7 +1302,7 @@ class CMFDRun(object): tally.filters = [mesh_filter, energy_filter] else: tally.filters = [mesh_filter] - # Set scores for tally + # Set scores, type, and estimator for tally tally.scores = ['flux', 'total'] tally.type = 'volume' tally.estimator = 'analog' @@ -1095,7 +1314,7 @@ class CMFDRun(object): tally.filters = [mesh_filter, energy_filter, energyout_filter] else: tally.filters = [mesh_filter] - # Set scores for tally + # Set scores, type, and estimator for tally tally.scores = ['nu-scatter', 'nu-fission'] tally.type = 'volume' tally.estimator = 'analog' @@ -1107,7 +1326,7 @@ class CMFDRun(object): tally.filters = [meshsurface_filter, energy_filter] else: tally.filters = [meshsurface_filter] - # Set scores for tally + # Set scores, type, and estimator for tally tally.scores = ['current'] tally.type = 'mesh-surface' tally.estimator = 'analog' diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ce4825426..38ad01428 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -254,7 +254,7 @@ contains ! Set the energy bin if needed if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + filter_matches(i_filter_ein) % bins % data(1) = 12*(ng - h) + 1 end if score_index = 0 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 60bed6425..d2d8eef02 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -65,7 +65,7 @@ module simulation_header ! ============================================================================ ! MISCELLANEOUS VARIABLES - integer :: restart_batch + integer(C_INT), bind(C, name='openmc_restart_batch') :: restart_batch ! Flag for enabling cell overlap checking during transport integer(8), allocatable :: overlap_check_cnt(:)