diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 3d2e8f57bf..0acc29e780 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -138,7 +138,7 @@ extern "C" { const int* indices, int n_elements, int dim, double spectral, const int* cmfd_indices, - const int* map); + const int* map, int n_threads); //! Runs a Gauss Seidel linear solver to solve CMFD matrix equations //! linear solver diff --git a/openmc/capi/core.py b/openmc/capi/core.py index a470f0665b..fe911d3a52 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -47,7 +47,7 @@ _dll.openmc_get_keff.argtypes = [POINTER(c_double*2)] _dll.openmc_get_keff.restype = c_int _dll.openmc_get_keff.errcheck = _error_handler _init_linsolver_argtypes = [_array_1d_int, c_int, _array_1d_int, c_int, c_int, - c_double, _array_1d_int, _array_1d_int] + c_double, _array_1d_int, _array_1d_int, c_int] _dll.openmc_initialize_linsolver.argtypes = _init_linsolver_argtypes _dll.openmc_initialize_linsolver.restype = None _dll.openmc_is_statepoint_batch.restype = c_bool diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 69a24ba617..c711c9b468 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -196,7 +196,7 @@ class CMFDRun(object): ---------- tally_begin : int Batch number at which CMFD tallies should begin accummulating - cmfd_begin: int + solver_begin: int Batch number at which CMFD solver should start executing ref_d : list of floats List of reference diffusion coefficients to fix CMFD parameters to @@ -307,8 +307,8 @@ class CMFDRun(object): """ # Variables that users can modify self._tally_begin = 1 - self._cmfd_begin = 1 - self._ref_d = [] + self._solver_begin = 1 + self._ref_d = np.array([]) self._display = {'balance': False, 'dominance': False, 'entropy': False, 'source': False} self._downscatter = False @@ -328,6 +328,7 @@ class CMFDRun(object): self._window_type = 'none' self._window_size = 10 self._intracomm = None + self._n_threads = 1 # External variables used during runtime but users cannot control self._set_reference_params = False @@ -413,8 +414,8 @@ class CMFDRun(object): return self._tally_begin @property - def cmfd_begin(self): - return self._cmfd_begin + def solver_begin(self): + return self._solver_begin @property def ref_d(self): @@ -492,6 +493,10 @@ class CMFDRun(object): def indices(self): return self._indices + @property + def n_threads(self): + return self._n_threads + @property def cmfd_src(self): return self._cmfd_src @@ -522,11 +527,12 @@ class CMFDRun(object): check_greater_than('CMFD tally begin batch', begin, 0) self._tally_begin = begin - @cmfd_begin.setter - def cmfd_begin(self, begin): + @solver_begin.setter + def solver_begin(self, begin): check_type('CMFD feedback begin batch', begin, Integral) check_greater_than('CMFD feedback begin batch', begin, 0) - self._cmfd_begin = begin + self._solver_begin = begin + @ref_d.setter def ref_d(self, diff_params): @@ -677,6 +683,12 @@ class CMFDRun(object): check_length('Gauss-Seidel tolerance', gauss_seidel_tolerance, 2) self._gauss_seidel_tolerance = gauss_seidel_tolerance + @n_threads.setter + def n_threads(self, n_threads): + check_type('CMFD number of threads', n_threads, Integral) + check_greater_than('CMFD number of threads', n_threads, 0) + self._n_threads = n_threads + def run(self, **kwargs): """Run OpenMC with coarse mesh finite difference acceleration @@ -692,7 +704,8 @@ class CMFDRun(object): """ with self.run_in_memory(**kwargs): for _ in self.iter_batches(): - pass + print('done') + #pass @contextmanager def run_in_memory(self, **kwargs): @@ -787,6 +800,7 @@ class CMFDRun(object): # Run next batch status = openmc.capi.next_batch() + print('2') # Perform CMFD calculations self._execute_cmfd() @@ -848,7 +862,7 @@ class CMFDRun(object): cmfd_group = f.create_group("cmfd") cmfd_group.attrs['cmfd_on'] = self._cmfd_on cmfd_group.attrs['feedback'] = self._feedback - cmfd_group.attrs['cmfd_begin'] = self._cmfd_begin + cmfd_group.attrs['solver_begin'] = self._solver_begin cmfd_group.attrs['mesh_id'] = self._mesh_id cmfd_group.attrs['tally_begin'] = self._tally_begin cmfd_group.attrs['time_cmfd'] = self._time_cmfd @@ -904,7 +918,7 @@ class CMFDRun(object): args = temp_loss.indptr, len(temp_loss.indptr), \ temp_loss.indices, len(temp_loss.indices), n, \ - self._spectral, self._indices, coremap + self._spectral, self._indices, coremap, self._n_threads return openmc.capi._dll.openmc_initialize_linsolver(*args) def _write_cmfd_output(self): @@ -933,9 +947,9 @@ class CMFDRun(object): """Write CMFD timing stats to buffer after finalizing simulation""" outstr = ("=====================> " "CMFD TIMING STATISTICS <====================\n\n" - " Time in CMFD = {:.5E} seconds\n" - " Building matrices = {:.5E} seconds\n" - " Solving matrices = {:.5E} seconds\n") + " Time in CMFD = {:.5e} seconds\n" + " Building matrices = {:.5e} seconds\n" + " Solving matrices = {:.5e} seconds\n") print(outstr.format(self._time_cmfd, self._time_cmfdbuild, self._time_cmfdsolve)) sys.stdout.flush() @@ -999,7 +1013,7 @@ class CMFDRun(object): dtype=int) # Check CMFD tallies accummulated before feedback turned on - if self._feedback and self._cmfd_begin < self._tally_begin: + if self._feedback and self._solver_begin < self._tally_begin: raise ValueError('Tally begin must be less than or equal to ' 'CMFD begin') @@ -1058,7 +1072,7 @@ class CMFDRun(object): # Define variables that exist on all processes self._cmfd_on = cmfd_group.attrs['cmfd_on'] self._feedback = cmfd_group.attrs['feedback'] - self._cmfd_begin = cmfd_group.attrs['cmfd_begin'] + self._solver_begin = cmfd_group.attrs['solver_begin'] self._tally_begin = cmfd_group.attrs['tally_begin'] self._k_cmfd = list(cmfd_group['k_cmfd']) self._dom = list(cmfd_group['dom']) @@ -1122,7 +1136,7 @@ class CMFDRun(object): current_batch = openmc.capi.current_batch() + 1 # Check to activate CMFD solver and possible feedback - if self._cmfd_begin == current_batch: + if self._solver_begin == current_batch: self._cmfd_on = True # Check to reset tallies @@ -1436,7 +1450,7 @@ class CMFDRun(object): source_energies = openmc.capi.source_bank()['E'] # Convert xyz location to the CMFD mesh index - mesh_ijk = np.floor((source_xyz-m.lower_left)/m.width).astype(int) + mesh_ijk = np.floor((source_xyz - m.lower_left)/m.width).astype(int) # Determine which energy bin each particle's energy belongs to # Separate into cases bases on where source energies lies on egrid @@ -1455,10 +1469,10 @@ class CMFDRun(object): mesh_ijk[:,0], mesh_ijk[:,1], mesh_ijk[:,2], energy_bins] if openmc.capi.master() and np.any(source_energies < energy[0]): - print(' WARNING: Source pt below energy grid') + print(' WARNING: Source point below energy grid') sys.stdout.flush() if openmc.capi.master() and np.any(source_energies > energy[-1]): - print(' WARNING: Source pt above energy grid') + print(' WARNING: Source point above energy grid') sys.stdout.flush() def _count_bank_sites(self): @@ -1814,8 +1828,8 @@ class CMFDRun(object): if self._power_monitor and openmc.capi.master(): str1 = ' {:d}:'.format(iter) str2 = 'k-eff: {:0.8f}'.format(k_n) - str3 = 'k-error: {:.5E}'.format(kerr) - str4 = 'src-error: {:.5E}'.format(serr) + str3 = 'k-error: {:.5e}'.format(kerr) + str4 = 'src-error: {:.5e}'.format(serr) str5 = ' {:d}'.format(innerits) print('{:8s}{:20s}{:25s}{:s}{:s}'.format(str1, str2, str3, str4, str5)) @@ -1898,11 +1912,12 @@ class CMFDRun(object): self._flux_rate = np.append(self._flux_rate, reshape_flux, axis=4) # Compute flux as aggregate of banked flux_rate over tally window - self._flux = np.where(is_accel[...,np.newaxis], + self._flux = np.where(is_accel[..., np.newaxis], np.sum(self._flux_rate, axis=4), 0.0) # Detect zero flux, abort if located and cmfd is on - zero_flux = np.logical_and(self._flux < _TINY_BIT, is_accel[...,np.newaxis]) + zero_flux = np.logical_and(self._flux < _TINY_BIT, + is_accel[..., np.newaxis]) if np.any(zero_flux) and self._cmfd_on: # Get index of first zero flux in flux array idx = np.argwhere(zero_flux)[0] @@ -2033,7 +2048,7 @@ class CMFDRun(object): axis=5) # Compute current as aggregate of banked current_rate over tally window - self._current = np.where(is_accel[...,np.newaxis,np.newaxis], + self._current = np.where(is_accel[..., np.newaxis, np.newaxis], np.sum(self._current_rate, axis=5), 0.0) # Get p1 scatter rr from CMFD tally 3 @@ -2142,12 +2157,12 @@ class CMFDRun(object): # Compute scattering rr by broadcasting flux in outgoing energy and # summing over incoming energy - scattering = np.sum(self._scattxs * self._flux[:,:,:,:,np.newaxis], + scattering = np.sum(self._scattxs * self._flux[:,:,:,:, np.newaxis], axis=3) # Compute fission rr by broadcasting flux in outgoing energy and # summing over incoming energy - fission = np.sum(self._nfissxs * self._flux[:,:,:,:,np.newaxis], + fission = np.sum(self._nfissxs * self._flux[:,:,:,:, np.newaxis], axis=3) # Compute residual @@ -2187,11 +2202,11 @@ class CMFDRun(object): self._dhat = np.zeros((nx, ny, nz, ng, 6)) # Set reference diffusion parameters - if list(self._ref_d): + if self._ref_d.size > 0: self._set_reference_params = True # Check length of reference diffusion parameters equal to number of # energy groups - if len(self._ref_d) != self._indices[3]: + if self._ref_d.size != self._indices[3]: raise OpenMCError('Number of reference diffusion parameters ' 'must equal number of CMFD energy groups') diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 27d594390c..40be440d5f 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -3,6 +3,9 @@ #include #include +#ifdef _OPENMP +#include +#endif #include "xtensor/xtensor.hpp" #include "openmc/error.h" @@ -29,6 +32,10 @@ int nx, ny, nz, ng; xt::xtensor indexmap; +int n_threads; + +int n_threads_reset; + } // namespace cmfd //============================================================================== @@ -101,7 +108,7 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x, for (int irb = 0; irb < 2; irb++) { // Loop around matrix rows - #pragma omp parallel for reduction (+:err) + #pragma omp parallel for reduction (+:err) num_threads(cmfd::n_threads) for (int irow = 0; irow < cmfd::dim; irow++) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -168,7 +175,7 @@ int cmfd_linsolver_2g(const double* A_data, const double* b, double* x, for (int irb = 0; irb < 2; irb++) { // Loop around matrix rows - #pragma omp parallel for reduction (+:err) + #pragma omp parallel for reduction (+:err) num_threads(cmfd::n_threads) for (int irow = 0; irow < cmfd::dim; irow+=2) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -304,7 +311,7 @@ extern "C" void openmc_initialize_linsolver(const int* indptr, int len_indptr, const int* indices, int n_elements, int dim, double spectral, const int* cmfd_indices, - const int* map) + const int* map, int n_threads) { // Store elements of indptr for (int i = 0; i < len_indptr; i++) @@ -331,6 +338,13 @@ void openmc_initialize_linsolver(const int* indptr, int len_indptr, cmfd::indexmap.resize({static_cast(dim), 3}); set_indexmap(map); } + +#ifdef _OPENMP + // Set number of threads to run CMFD solver on and store number of threads + // to reset to after solver finishes executing + cmfd::n_threads = n_threads; + cmfd::n_threads_reset = omp_get_max_threads(); +#endif } //============================================================================== @@ -342,14 +356,20 @@ extern "C" int openmc_run_linsolver(const double* A_data, const double* b, double* x, double tol) { + int result; + switch (cmfd::ng) { case 1: - return cmfd_linsolver_1g(A_data, b, x, tol); + result = cmfd_linsolver_1g(A_data, b, x, tol); + break; case 2: - return cmfd_linsolver_2g(A_data, b, x, tol); + result = cmfd_linsolver_2g(A_data, b, x, tol); + break; default: - return cmfd_linsolver_ng(A_data, b, x, tol); + result = cmfd_linsolver_ng(A_data, b, x, tol); + break; } + return result; } void free_memory_cmfd() diff --git a/tests/regression_tests/cmfd_feed/test.py b/tests/regression_tests/cmfd_feed/test.py index 87370a894f..6b56635f4a 100644 --- a/tests/regression_tests/cmfd_feed/test.py +++ b/tests/regression_tests/cmfd_feed/test.py @@ -24,7 +24,7 @@ def test_cmfd_physical_adjoint(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.run_adjoint = True @@ -54,7 +54,7 @@ def test_cmfd_math_adjoint(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.run_adjoint = True @@ -83,7 +83,7 @@ def test_cmfd_write_matrices(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.display = {'dominance': True} cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] @@ -131,7 +131,7 @@ def test_cmfd_feed(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.display = {'dominance': True} cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] diff --git a/tests/regression_tests/cmfd_feed_2g/test.py b/tests/regression_tests/cmfd_feed_2g/test.py index 8f0fcda686..d3af8998b6 100644 --- a/tests/regression_tests/cmfd_feed_2g/test.py +++ b/tests/regression_tests/cmfd_feed_2g/test.py @@ -17,7 +17,7 @@ def test_cmfd_feed_2g(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.display = {'dominance': True} cmfd_run.feedback = True cmfd_run.downscatter = True diff --git a/tests/regression_tests/cmfd_feed_expanding_window/test.py b/tests/regression_tests/cmfd_feed_expanding_window/test.py index 35cf3c2898..964d4f2253 100644 --- a/tests/regression_tests/cmfd_feed_expanding_window/test.py +++ b/tests/regression_tests/cmfd_feed_expanding_window/test.py @@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 10 + cmfd_run.solver_begin = 10 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.window_type = 'expanding' diff --git a/tests/regression_tests/cmfd_feed_ng/test.py b/tests/regression_tests/cmfd_feed_ng/test.py index 6b9ff60c56..a2a522e9c4 100644 --- a/tests/regression_tests/cmfd_feed_ng/test.py +++ b/tests/regression_tests/cmfd_feed_ng/test.py @@ -18,7 +18,7 @@ def test_cmfd_feed_ng(): cmfd_run.mesh = cmfd_mesh cmfd_run.reset = [5] cmfd_run.tally_begin = 10 - cmfd_run.cmfd_begin = 10 + cmfd_run.solver_begin = 10 cmfd_run.display = {'dominance': True} cmfd_run.feedback = True cmfd_run.downscatter = True diff --git a/tests/regression_tests/cmfd_feed_ref_d/test.py b/tests/regression_tests/cmfd_feed_ref_d/test.py index 925fec32d5..120d94b6b4 100644 --- a/tests/regression_tests/cmfd_feed_ref_d/test.py +++ b/tests/regression_tests/cmfd_feed_ref_d/test.py @@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 10 + cmfd_run.solver_begin = 10 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.window_type = 'expanding' diff --git a/tests/regression_tests/cmfd_feed_rolling_window/test.py b/tests/regression_tests/cmfd_feed_rolling_window/test.py index 07eb0741bb..2c7b7f242c 100644 --- a/tests/regression_tests/cmfd_feed_rolling_window/test.py +++ b/tests/regression_tests/cmfd_feed_rolling_window/test.py @@ -15,7 +15,7 @@ def test_cmfd_feed_rolling_window(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 10 + cmfd_run.solver_begin = 10 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.window_type = 'rolling' diff --git a/tests/regression_tests/cmfd_nofeed/test.py b/tests/regression_tests/cmfd_nofeed/test.py index 35c18a38f9..7ab72f86ad 100644 --- a/tests/regression_tests/cmfd_nofeed/test.py +++ b/tests/regression_tests/cmfd_nofeed/test.py @@ -15,7 +15,7 @@ def test_cmfd_nofeed(): # Initialize and run CMFDRun object cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.display = {'dominance': True} cmfd_run.feedback = False cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] diff --git a/tests/regression_tests/cmfd_restart/test.py b/tests/regression_tests/cmfd_restart/test.py index 369ef0349f..8f8410a7aa 100644 --- a/tests/regression_tests/cmfd_restart/test.py +++ b/tests/regression_tests/cmfd_restart/test.py @@ -53,7 +53,7 @@ def test_cmfd_restart(): cmfd_run = cmfd.CMFDRun() cmfd_run.mesh = cmfd_mesh cmfd_run.tally_begin = 5 - cmfd_run.cmfd_begin = 5 + cmfd_run.solver_begin = 5 cmfd_run.feedback = True cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20] cmfd_run.run() @@ -62,7 +62,7 @@ def test_cmfd_restart(): cmfd_run2 = cmfd.CMFDRun() cmfd_run2.mesh = cmfd_mesh2 cmfd_run2.tally_begin = 5 - cmfd_run2.cmfd_begin = 5 + cmfd_run2.solver_begin = 5 cmfd_run2.feedback = True cmfd_run2.gauss_seidel_tolerance = [1.e-15, 1.e-20]