From 0da3bf338180255d80a3a997bebe8445d4667a82 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 9 Jul 2019 19:54:40 -0400 Subject: [PATCH 01/16] Try OMP parallelism --- src/cmfd_solver.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index f7ef57546f..c46ace3f73 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -101,6 +101,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 for (int irow = 0; irow < cmfd::dim; irow++) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -167,6 +168,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 for (int irow = 0; irow < cmfd::dim; irow+=2) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -255,6 +257,7 @@ int cmfd_linsolver_ng(const double* A_data, const double* b, double* x, std::vector tmpx {x, x+cmfd::dim}; // Loop around matrix rows +#pragma omp parallel for for (int irow = 0; irow < cmfd::dim; irow++) { // Get index of diagonal for current row int didx = get_diagonal_index(irow); From 3974631b7255d87416ba83cccd99823c69e6e9a8 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Wed, 10 Jul 2019 13:00:29 -0400 Subject: [PATCH 02/16] reduce err across threads, modify xsdata scripts --- scripts/openmc-make-test-data | 2 +- src/cmfd_solver.cpp | 10 +++++----- tools/ci/download-xs.sh | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/scripts/openmc-make-test-data b/scripts/openmc-make-test-data index a37a78fa9e..439dcde6a3 100755 --- a/scripts/openmc-make-test-data +++ b/scripts/openmc-make-test-data @@ -158,7 +158,7 @@ with tempfile.TemporaryDirectory() as tmpdir: print('Creating compressed archive...') test_tar = pwd / 'nndc_hdf5_test.tar.xz' with tarfile.open(str(test_tar), 'w:xz') as txz: - txz.add('nndc_hdf5') + txz.add(output_dir) # Change back to original directory os.chdir(str(pwd)) diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index c46ace3f73..472107d0cb 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -101,7 +101,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 + #pragma omp parallel for reduction (+:err) for (int irow = 0; irow < cmfd::dim; irow++) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -127,12 +127,13 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x, // Compute residual and update error double res = (tmpx[irow] - x[irow]) / tmpx[irow]; - err += res * res; + err = res * res; } } // Check convergence err = std::sqrt(err / cmfd::dim); + std::cout << err << "\n"; if (err < tol) return igs; @@ -168,7 +169,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 + #pragma omp parallel for reduction (+:err) for (int irow = 0; irow < cmfd::dim; irow+=2) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -219,7 +220,7 @@ int cmfd_linsolver_2g(const double* A_data, const double* b, double* x, // Compute residual and update error double res = (tmpx[irow] - x[irow]) / tmpx[irow]; - err += res * res; + err = res * res; } } @@ -257,7 +258,6 @@ int cmfd_linsolver_ng(const double* A_data, const double* b, double* x, std::vector tmpx {x, x+cmfd::dim}; // Loop around matrix rows -#pragma omp parallel for for (int irow = 0; irow < cmfd::dim; irow++) { // Get index of diagonal for current row int didx = get_diagonal_index(irow); diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 07ade9c70f..e831d8215a 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -7,7 +7,7 @@ if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then fi # Download ENDF/B-VII.1 distribution -ENDF=$HOME/endf-b-vii.1/ +ENDF=$HOME/endf-b-vii.1 if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ fi From 596eab4bb69743a85e1f03830229eedb1091957d Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 22 Jul 2019 12:08:04 -0400 Subject: [PATCH 03/16] Define default value for resnb when flux equals 0 --- openmc/cmfd.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 742cdeabfd..ca52502c3a 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1141,7 +1141,7 @@ class CMFDRun(object): self._dhat = np.zeros((nx, ny, nz, ng, 6)) # Set reference diffusion parameters - if self._ref_d: + if list(self._ref_d): self._set_reference_params = True # Check length of reference diffusion parameters equal to number of # energy groups @@ -2214,7 +2214,8 @@ class CMFDRun(object): res = leakage + interactions - scattering - (1.0 / keff) * fission # Normalize res by flux and bank res - self._resnb = np.divide(res, self._flux, where=self._flux > 0) + self._resnb = np.divide(res, self._flux, where=self._flux > 0, + out=np.zeros_like(self._flux)) # Calculate RMS and record for this batch self._balance.append(np.sqrt( From 6e77ddc532f4bfb322c57f93bd58411d4c248f64 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 22 Jul 2019 14:07:58 -0400 Subject: [PATCH 04/16] Reduce err properly --- src/cmfd_solver.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 472107d0cb..11d6d97b5d 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -133,7 +133,6 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x, // Check convergence err = std::sqrt(err / cmfd::dim); - std::cout << err << "\n"; if (err < tol) return igs; @@ -220,7 +219,7 @@ int cmfd_linsolver_2g(const double* A_data, const double* b, double* x, // Compute residual and update error double res = (tmpx[irow] - x[irow]) / tmpx[irow]; - err = res * res; + err += res * res; } } From bc363cc934453e423901462c43616a4b74558b22 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 22 Jul 2019 14:14:49 -0400 Subject: [PATCH 05/16] Apply changes to 1g solver --- src/cmfd_solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 11d6d97b5d..27d594390c 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -127,7 +127,7 @@ int cmfd_linsolver_1g(const double* A_data, const double* b, double* x, // Compute residual and update error double res = (tmpx[irow] - x[irow]) / tmpx[irow]; - err = res * res; + err += res * res; } } From 093f90fc6460fd2be7f0e590ab0d076cfb4cce82 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 22 Jul 2019 16:51:29 -0400 Subject: [PATCH 06/16] Change logic to run CMFD at cmfd_begin instead of tally_begin --- openmc/cmfd.py | 223 +++++++++--------- tests/regression_tests/cmfd_feed/test.py | 8 +- tests/regression_tests/cmfd_feed_2g/test.py | 2 +- .../results_true.dat | 25 -- .../cmfd_feed_expanding_window/test.py | 2 +- tests/regression_tests/cmfd_feed_ng/test.py | 2 +- .../cmfd_feed_ref_d/results_true.dat | 25 -- .../regression_tests/cmfd_feed_ref_d/test.py | 2 +- .../cmfd_feed_rolling_window/results_true.dat | 25 -- .../cmfd_feed_rolling_window/test.py | 2 +- tests/regression_tests/cmfd_nofeed/test.py | 2 +- tests/regression_tests/cmfd_restart/test.py | 4 +- 12 files changed, 120 insertions(+), 202 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index ca52502c3a..8fa406bc5a 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -196,8 +196,8 @@ class CMFDRun(object): ---------- tally_begin : int Batch number at which CMFD tallies should begin accummulating - feedback_begin: int - Batch number at which CMFD feedback should be turned on + cmfd_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 dhat_reset : bool @@ -310,7 +310,7 @@ class CMFDRun(object): """ # Variables that users can modify self._tally_begin = 1 - self._feedback_begin = 1 + self._cmfd_begin = 1 self._ref_d = [] self._dhat_reset = False self._display = {'balance': False, 'dominance': False, @@ -418,8 +418,8 @@ class CMFDRun(object): return self._tally_begin @property - def feedback_begin(self): - return self._feedback_begin + def cmfd_begin(self): + return self._cmfd_begin @property def ref_d(self): @@ -531,11 +531,11 @@ class CMFDRun(object): check_greater_than('CMFD tally begin batch', begin, 0) self._tally_begin = begin - @feedback_begin.setter - def feedback_begin(self, begin): + @cmfd_begin.setter + def cmfd_begin(self, begin): check_type('CMFD feedback begin batch', begin, Integral) check_greater_than('CMFD feedback begin batch', begin, 0) - self._feedback_begin = begin + self._cmfd_begin = begin @ref_d.setter def ref_d(self, diff_params): @@ -801,13 +801,8 @@ class CMFDRun(object): # Run next batch status = openmc.capi.next_batch() - # Perform CMFD calculation if on - if self._cmfd_on: - self._execute_cmfd() - - # Write CMFD output if CMFD on for current batch - if openmc.capi.master(): - self._write_cmfd_output() + # Perform CMFD calculations + self._execute_cmfd() # Write CMFD data to statepoint if openmc.capi.is_statepoint_batch(): @@ -865,7 +860,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['feedback_begin'] = self._feedback_begin + cmfd_group.attrs['cmfd_begin'] = self._cmfd_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 @@ -1040,9 +1035,9 @@ class CMFDRun(object): dtype=int) # Check CMFD tallies accummulated before feedback turned on - if self._feedback and self._feedback_begin < self._tally_begin: + if self._feedback and self._cmfd_begin < self._tally_begin: raise ValueError('Tally begin must be less than or equal to ' - 'feedback begin') + 'CMFD begin') # Set number of batches where cmfd tallies should be reset self._n_resets = len(self._reset) @@ -1072,7 +1067,7 @@ class CMFDRun(object): cmfd_group = f['cmfd'] self._cmfd_on = cmfd_group.attrs['cmfd_on'] self._feedback = cmfd_group.attrs['feedback'] - self._feedback_begin = cmfd_group.attrs['feedback_begin'] + self._cmfd_begin = cmfd_group.attrs['cmfd_begin'] self._tally_begin = cmfd_group.attrs['tally_begin'] self._time_cmfd = cmfd_group.attrs['time_cmfd'] self._time_cmfdbuild = cmfd_group.attrs['time_cmfdbuild'] @@ -1169,9 +1164,8 @@ class CMFDRun(object): # Add 1 as next_batch has not been called yet current_batch = openmc.capi.current_batch() + 1 - # Check to activate CMFD diffusion and possible feedback - # Check to activate CMFD tallies - if self._tally_begin == current_batch: + # Check to activate CMFD solver and possible feedback + if self._cmfd_begin == current_batch: self._cmfd_on = True # Check to reset tallies @@ -1181,35 +1175,46 @@ class CMFDRun(object): def _execute_cmfd(self): """Runs CMFD calculation on master node""" - # Run CMFD on single processor on master if openmc.capi.master(): # Start CMFD timer time_start_cmfd = time.time() - # Create CMFD data from OpenMC tallies - self._set_up_cmfd() + if openmc.capi.current_batch() >= self._tally_begin: + # Calculate all cross sections based on tally window averages + self._compute_xs() - # Call solver - self._cmfd_solver_execute() + # Execute CMFD algorithm if CMFD on for current batch + if self._cmfd_on: + # Run CMFD on single processor on master + if openmc.capi.master(): + # Create CMFD data based on OpenMC tallies + self._set_up_cmfd() - # Store k-effective - self._k_cmfd.append(self._keff) + # Call solver + self._cmfd_solver_execute() - # Check to perform adjoint on last batch - if (openmc.capi.current_batch() == openmc.capi.settings.batches - and self._run_adjoint): - self._cmfd_solver_execute(adjoint=True) + # Store k-effective + self._k_cmfd.append(self._keff) - # Calculate fission source - self._calc_fission_source() + # Check to perform adjoint on last batch + if (openmc.capi.current_batch() == openmc.capi.settings.batches + and self._run_adjoint): + self._cmfd_solver_execute(adjoint=True) - # Calculate weight factors - self._cmfd_reweight(True) + # Calculate fission source + self._calc_fission_source() + + # Calculate weight factors + self._cmfd_reweight() # Stop CMFD timer if openmc.capi.master(): time_stop_cmfd = time.time() self._time_cmfd += time_stop_cmfd - time_start_cmfd + if self._cmfd_on: + # Write CMFD output if CMFD on for current batch + self._write_cmfd_output() + def _cmfd_tally_reset(self): """Resets all CMFD tallies in memory""" @@ -1228,9 +1233,6 @@ class CMFDRun(object): """Configures CMFD object for a CMFD eigenvalue calculation """ - # Calculate all cross sections based on tally window averages - self._compute_xs() - # Compute effective downscatter cross section if self._downscatter: self._compute_effective_downscatter() @@ -1422,96 +1424,85 @@ class CMFDRun(object): self._src_cmp.append(np.sqrt(1.0 / self._norm * np.sum((self._cmfd_src - self._openmc_src)**2))) - def _cmfd_reweight(self, new_weights): - """Performs weighting of particles in source bank + def _cmfd_reweight(self): + """Performs weighting of particles in source bank""" + # Get spatial dimensions and energy groups + nx, ny, nz, ng = self._indices - Parameters - ---------- - new_weights : bool - Whether to reweight particles or not + # Count bank site in mesh and reverse due to egrid structured + outside = self._count_bank_sites() - """ - # Compute new weight factors - if new_weights: + # Check and raise error if source sites exist outside of CMFD mesh + if openmc.capi.master() and outside: + raise OpenMCError('Source sites outside of the CMFD mesh') - # Get spatial dimensions and energy groups - nx, ny, nz, ng = self._indices + # Have master compute weight factors, ignore any zeros in + # sourcecounts or cmfd_src + if openmc.capi.master(): + # Compute normalization factor + norm = np.sum(self._sourcecounts) / np.sum(self._cmfd_src) - # Count bank site in mesh and reverse due to egrid structured - outside = self._count_bank_sites() + # Define target reshape dimensions for sourcecounts. This + # defines how self._sourcecounts is ordered by dimension + target_shape = [nz, ny, nx, ng] - # Check and raise error if source sites exist outside of CMFD mesh - if openmc.capi.master() and outside: - raise OpenMCError('Source sites outside of the CMFD mesh') + # Reshape sourcecounts to target shape. Swap x and z axes so + # that the shape is now [nx, ny, nz, ng] + sourcecounts = np.swapaxes( + self._sourcecounts.reshape(target_shape), 0, 2) - # Have master compute weight factors, ignore any zeros in - # sourcecounts or cmfd_src - if openmc.capi.master(): - # Compute normalization factor - norm = np.sum(self._sourcecounts) / np.sum(self._cmfd_src) + # Flip index of energy dimension + sourcecounts = np.flip(sourcecounts, axis=3) - # Define target reshape dimensions for sourcecounts. This - # defines how self._sourcecounts is ordered by dimension - target_shape = [nz, ny, nx, ng] + # Compute weight factors + div_condition = np.logical_and(sourcecounts > 0, + self._cmfd_src > 0) + self._weightfactors = (np.divide(self._cmfd_src * norm, + sourcecounts, where=div_condition, + out=np.ones_like(self._cmfd_src), + dtype=np.float32)) - # Reshape sourcecounts to target shape. Swap x and z axes so - # that the shape is now [nx, ny, nz, ng] - sourcecounts = np.swapaxes( - self._sourcecounts.reshape(target_shape), 0, 2) + if not self._feedback: + return - # Flip index of energy dimension - sourcecounts = np.flip(sourcecounts, axis=3) + # Broadcast weight factors to all procs + if have_mpi: + self._weightfactors = self._intracomm.bcast( + self._weightfactors) - # Compute weight factors - div_condition = np.logical_and(sourcecounts > 0, - self._cmfd_src > 0) - self._weightfactors = (np.divide(self._cmfd_src * norm, - sourcecounts, where=div_condition, - out=np.ones_like(self._cmfd_src), - dtype=np.float32)) + m = openmc.capi.meshes[self._mesh_id] + energy = self._egrid + ng = self._indices[3] - if (not self._feedback - or openmc.capi.current_batch() < self._feedback_begin): - return + # Get locations and energies of all particles in source bank + source_xyz = openmc.capi.source_bank()['r'] + source_energies = openmc.capi.source_bank()['E'] - # Broadcast weight factors to all procs - if have_mpi: - self._weightfactors = self._intracomm.bcast( - self._weightfactors) + # Convert xyz location to the CMFD mesh index + mesh_ijk = np.floor((source_xyz-m.lower_left)/m.width).astype(int) - m = openmc.capi.meshes[self._mesh_id] - energy = self._egrid - ng = self._indices[3] + # Determine which energy bin each particle's energy belongs to + # Separate into cases bases on where source energies lies on egrid + energy_bins = np.zeros(len(source_energies), dtype=int) + idx = np.where(source_energies < energy[0]) + energy_bins[idx] = ng - 1 + idx = np.where(source_energies > energy[-1]) + energy_bins[idx] = 0 + idx = np.where((source_energies >= energy[0]) & + (source_energies <= energy[-1])) + energy_bins[idx] = ng - np.digitize(source_energies, energy) - # Get locations and energies of all particles in source bank - source_xyz = openmc.capi.source_bank()['r'] - source_energies = openmc.capi.source_bank()['E'] + # Determine weight factor of each particle based on its mesh index + # and energy bin and updates its weight + openmc.capi.source_bank()['wgt'] *= self._weightfactors[ + mesh_ijk[:,0], mesh_ijk[:,1], mesh_ijk[:,2], energy_bins] - # Convert xyz location to the CMFD mesh index - 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 - energy_bins = np.zeros(len(source_energies), dtype=int) - idx = np.where(source_energies < energy[0]) - energy_bins[idx] = ng - 1 - idx = np.where(source_energies > energy[-1]) - energy_bins[idx] = 0 - idx = np.where((source_energies >= energy[0]) & - (source_energies <= energy[-1])) - energy_bins[idx] = ng - np.digitize(source_energies, energy) - - # Determine weight factor of each particle based on its mesh index - # and energy bin and updates its weight - openmc.capi.source_bank()['wgt'] *= self._weightfactors[ - 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') - sys.stdout.flush() - if openmc.capi.master() and np.any(source_energies > energy[-1]): - print(' WARNING: Source pt above energy grid') - sys.stdout.flush() + if openmc.capi.master() and np.any(source_energies < energy[0]): + print(' WARNING: Source pt below energy grid') + sys.stdout.flush() + if openmc.capi.master() and np.any(source_energies > energy[-1]): + print(' WARNING: Source pt above energy grid') + sys.stdout.flush() def _count_bank_sites(self): """Determines the number of fission bank sites in each cell of a given @@ -1939,6 +1930,8 @@ class CMFDRun(object): tally_results = tallies[tally_id].results[:,0,1] flux = np.where(is_cmfd_accel, tally_results, 0.) + # TODO do this check after flux reshape + # TODO need to update is_cmfd_accel, current, and coremap # Detect zero flux, abort if located if np.any(flux[is_cmfd_accel] < _TINY_BIT): # Get index of zero flux in flux array diff --git a/tests/regression_tests/cmfd_feed/test.py b/tests/regression_tests/cmfd_feed/test.py index 906c631ecc..87370a894f 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.feedback_begin = 5 + cmfd_run.cmfd_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.feedback_begin = 5 + cmfd_run.cmfd_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.feedback_begin = 5 + cmfd_run.cmfd_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.feedback_begin = 5 + cmfd_run.cmfd_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 ba4d21609a..8f0fcda686 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.feedback_begin = 5 + cmfd_run.cmfd_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/results_true.dat b/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat index b67bef0f76..98b7a5e823 100644 --- a/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat +++ b/tests/regression_tests/cmfd_feed_expanding_window/results_true.dat @@ -391,11 +391,6 @@ cmfd indices 1.000000E+00 1.000000E+00 k cmfd -1.143597E+00 -1.163387E+00 -1.173384E+00 -1.171035E+00 -1.147196E+00 1.122260E+00 1.106380E+00 1.124693E+00 @@ -408,11 +403,6 @@ k cmfd 1.170308E+00 1.184540E+00 cmfd entropy -3.212002E+00 -3.206393E+00 -3.223984E+00 -3.222764E+00 -3.232123E+00 3.242083E+00 3.246067E+00 3.238869E+00 @@ -425,11 +415,6 @@ cmfd entropy 3.221485E+00 3.219108E+00 cmfd balance -8.26180E-03 -4.27338E-03 -2.22686E-03 -1.93026E-03 -1.96979E-03 2.13756E-03 2.01479E-03 1.74519E-03 @@ -442,11 +427,6 @@ cmfd balance 1.24780E-03 1.15560E-03 cmfd dominance ratio -5.404E-01 -5.406E-01 -5.449E-01 -5.473E-01 -5.534E-01 5.623E-01 5.738E-01 5.611E-01 @@ -459,11 +439,6 @@ cmfd dominance ratio 5.412E-01 5.383E-01 cmfd openmc source comparison -1.575499E-02 -1.293688E-02 -3.531746E-03 -8.281178E-03 -5.771681E-03 7.459013E-03 5.012869E-03 1.770224E-03 diff --git a/tests/regression_tests/cmfd_feed_expanding_window/test.py b/tests/regression_tests/cmfd_feed_expanding_window/test.py index e8671e95c7..35cf3c2898 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.feedback_begin = 10 + cmfd_run.cmfd_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 dce674f905..6b9ff60c56 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.feedback_begin = 10 + cmfd_run.cmfd_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/results_true.dat b/tests/regression_tests/cmfd_feed_ref_d/results_true.dat index 06440d81f6..4a26ded630 100644 --- a/tests/regression_tests/cmfd_feed_ref_d/results_true.dat +++ b/tests/regression_tests/cmfd_feed_ref_d/results_true.dat @@ -391,11 +391,6 @@ cmfd indices 1.000000E+00 1.000000E+00 k cmfd -1.143785E+00 -1.163460E+00 -1.173453E+00 -1.171056E+00 -1.147214E+00 1.122230E+00 1.106385E+00 1.124706E+00 @@ -408,11 +403,6 @@ k cmfd 1.170183E+00 1.184408E+00 cmfd entropy -3.211758E+00 -3.206356E+00 -3.223933E+00 -3.222754E+00 -3.232110E+00 3.242098E+00 3.246062E+00 3.238858E+00 @@ -425,11 +415,6 @@ cmfd entropy 3.221498E+00 3.219126E+00 cmfd balance -8.26180E-03 -4.27338E-03 -2.22686E-03 -1.93026E-03 -1.96979E-03 2.13756E-03 2.01521E-03 1.74538E-03 @@ -442,11 +427,6 @@ cmfd balance 1.24170E-03 1.15645E-03 cmfd dominance ratio -5.408E-01 -5.374E-01 -5.420E-01 -5.444E-01 -5.513E-01 5.613E-01 5.722E-01 5.592E-01 @@ -459,11 +439,6 @@ cmfd dominance ratio 5.413E-01 5.381E-01 cmfd openmc source comparison -1.597982E-02 -1.276493E-02 -3.495229E-03 -8.157807E-03 -5.715330E-03 7.433111E-03 5.006211E-03 1.766072E-03 diff --git a/tests/regression_tests/cmfd_feed_ref_d/test.py b/tests/regression_tests/cmfd_feed_ref_d/test.py index d033ab8ddf..925fec32d5 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.feedback_begin = 10 + cmfd_run.cmfd_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/results_true.dat b/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat index 5180526835..98a8f8d06a 100644 --- a/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat +++ b/tests/regression_tests/cmfd_feed_rolling_window/results_true.dat @@ -391,11 +391,6 @@ cmfd indices 1.000000E+00 1.000000E+00 k cmfd -1.143597E+00 -1.163387E+00 -1.162391E+00 -1.163351E+00 -1.145721E+00 1.134785E+00 1.119048E+00 1.116124E+00 @@ -408,11 +403,6 @@ k cmfd 1.166709E+00 1.167223E+00 cmfd entropy -3.212002E+00 -3.206393E+00 -3.222109E+00 -3.221996E+00 -3.229978E+00 3.234615E+00 3.246512E+00 3.244634E+00 @@ -425,11 +415,6 @@ cmfd entropy 3.203758E+00 3.201798E+00 cmfd balance -8.26180E-03 -4.27338E-03 -2.62159E-03 -2.40301E-03 -2.08484E-03 1.58351E-03 1.59196E-03 1.87591E-03 @@ -442,11 +427,6 @@ cmfd balance 2.25515E-03 1.53613E-03 cmfd dominance ratio -5.404E-01 -5.406E-01 -5.432E-01 -5.454E-01 -5.507E-01 5.578E-01 5.679E-01 5.671E-01 @@ -459,11 +439,6 @@ cmfd dominance ratio 5.374E-01 5.321E-01 cmfd openmc source comparison -1.575499E-02 -1.293688E-02 -5.920734E-03 -9.746850E-03 -7.183896E-03 7.693485E-03 4.158805E-03 2.962505E-03 diff --git a/tests/regression_tests/cmfd_feed_rolling_window/test.py b/tests/regression_tests/cmfd_feed_rolling_window/test.py index 802b5b7003..07eb0741bb 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.feedback_begin = 10 + cmfd_run.cmfd_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 7d0895c2fb..35c18a38f9 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.tally_begin = 5 + cmfd_run.cmfd_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 7ed92d024c..369ef0349f 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.feedback_begin = 5 + cmfd_run.cmfd_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.feedback_begin = 5 + cmfd_run2.cmfd_begin = 5 cmfd_run2.feedback = True cmfd_run2.gauss_seidel_tolerance = [1.e-15, 1.e-20] From 72bdcbb48bdfdea1b6b03746a0a6b958da4d7421 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 22 Jul 2019 23:52:30 -0400 Subject: [PATCH 07/16] Thrown zero flux error only if sum over window is zero instead of at specific tally realization --- openmc/cmfd.py | 51 ++++++++++++++++++++------------------------------ 1 file changed, 20 insertions(+), 31 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 8fa406bc5a..ac425d2832 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1918,35 +1918,12 @@ class CMFDRun(object): # Get tallies in-memory tallies = openmc.capi.tallies - # Ravel coremap as 1d array similar to how tally data is arranged - coremap = np.ravel(self._coremap.swapaxes(0, 2)) - # Set conditional numpy array as boolean vector based on coremap - # Repeat each value for number of groups in problem - is_cmfd_accel = np.repeat(coremap != _CMFD_NOACCEL, ng) + is_accel = self._coremap != _CMFD_NOACCEL # Get flux from CMFD tally 0 tally_id = self._tally_ids[0] - tally_results = tallies[tally_id].results[:,0,1] - flux = np.where(is_cmfd_accel, tally_results, 0.) - - # TODO do this check after flux reshape - # TODO need to update is_cmfd_accel, current, and coremap - # 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 OpenMCError(err_message) + flux = tallies[tally_id].results[:,0,1] # Define target tally reshape dimensions. This defines how openmc # tallies are ordered by dimension @@ -1964,7 +1941,21 @@ 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.sum(self._flux_rate, axis=4) + 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 + if np.any(self._flux[is_accel[...,:]] < _TINY_BIT) and self.cmfd_on: + # Get index of first zero flux in flux array + idx = np.argwhere(self._flux[is_accel[...,:]] < _TINY_BIT)[0] + + # Throw error message (one-based indexing) + # Index of group is flipped + err_message = 'Detected zero flux without coremap overlay' + \ + ' at mesh: (' + \ + ', '.join(str(i+1) for i in idx[:-1]) + \ + ') in group ' + str(ng-idx[-1]) + raise OpenMCError(err_message) # Get total rr from CMFD tally 0 totalrr = tallies[tally_id].results[:,1,1] @@ -2065,10 +2056,7 @@ class CMFDRun(object): # Get surface currents from CMFD tally 2 tally_id = self._tally_ids[2] - tally_results = tallies[tally_id].results[:,0,1] - - # Filter tally results to include only accelerated regions - current = np.where(np.repeat(is_cmfd_accel, 12), tally_results, 0.) + current = tallies[tally_id].results[:,0,1] # Define target tally reshape dimensions for current target_tally_shape = [nz, ny, nx, 12, ng, 1] @@ -2087,7 +2075,8 @@ class CMFDRun(object): axis=5) # Compute current as aggregate of banked current_rate over tally window - self._current = np.sum(self._current_rate, axis=5) + self._current = np.where(is_accel[...,np.newaxis,np.newaxis], + np.sum(self._current_rate, axis=5), 0.0) # Get p1 scatter rr from CMFD tally 3 tally_id = self._tally_ids[3] From 11b3a1fe7be3f64578480f8e794f4dd99f614072 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 23 Jul 2019 11:27:22 -0400 Subject: [PATCH 08/16] Fix digitize bug --- openmc/cmfd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index ac425d2832..0a18dffdfb 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1490,7 +1490,7 @@ class CMFDRun(object): energy_bins[idx] = 0 idx = np.where((source_energies >= energy[0]) & (source_energies <= energy[-1])) - energy_bins[idx] = ng - np.digitize(source_energies, energy) + energy_bins[idx] = ng - np.digitize(source_energies[idx], energy) # Determine weight factor of each particle based on its mesh index # and energy bin and updates its weight @@ -1547,7 +1547,7 @@ class CMFDRun(object): energy_bins[idx] = ng - 1 idx = np.where((source_energies >= energy[0]) & (source_energies <= energy[-1])) - energy_bins[idx] = np.digitize(source_energies, energy) - 1 + energy_bins[idx] = np.digitize(source_energies[idx], energy) - 1 # Determine all unique combinations of mesh bin and energy bin, and # count number of particles that belong to these combinations From f3286c6601af99694f4c311c526f91b41d1ee1f1 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 23 Jul 2019 15:17:59 -0400 Subject: [PATCH 09/16] Fix typo to make cmfd_on private --- openmc/cmfd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index 0a18dffdfb..e4976c92c2 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1945,7 +1945,7 @@ class CMFDRun(object): np.sum(self._flux_rate, axis=4), 0.0) # Detect zero flux, abort if located and cmfd is on - if np.any(self._flux[is_accel[...,:]] < _TINY_BIT) and self.cmfd_on: + if np.any(self._flux[is_accel[...,:]] < _TINY_BIT) and self._cmfd_on: # Get index of first zero flux in flux array idx = np.argwhere(self._flux[is_accel[...,:]] < _TINY_BIT)[0] From 0f1b0a2ac6a12e3eff247815b104627747e4f2ac Mon Sep 17 00:00:00 2001 From: shikhark Date: Tue, 23 Jul 2019 21:18:16 +0000 Subject: [PATCH 10/16] Get proper indices for zero flux error --- openmc/cmfd.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index e4976c92c2..b0a4f55b18 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -1945,9 +1945,10 @@ class CMFDRun(object): np.sum(self._flux_rate, axis=4), 0.0) # Detect zero flux, abort if located and cmfd is on - if np.any(self._flux[is_accel[...,:]] < _TINY_BIT) and self._cmfd_on: + zero_flux = np.logical_and(self._flux < _TINY_BIT, is_accel[...,np.newaxis]) + if np.any(zero_flux) and self._cmfd_on: # Get index of first zero flux in flux array - idx = np.argwhere(self._flux[is_accel[...,:]] < _TINY_BIT)[0] + idx = np.argwhere(zero_flux)[0] # Throw error message (one-based indexing) # Index of group is flipped From c6adc5f0536483852462e41d940be0bad0b48c4c Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 30 Jul 2019 14:48:07 -0400 Subject: [PATCH 11/16] Separate variables defined on all processes vs just on master --- openmc/cmfd.py | 258 +++++++++++++++++++++++-------------------------- 1 file changed, 120 insertions(+), 138 deletions(-) diff --git a/openmc/cmfd.py b/openmc/cmfd.py index b0a4f55b18..69a24ba617 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -200,9 +200,6 @@ class CMFDRun(object): Batch number at which CMFD solver should start executing ref_d : list of floats List of reference diffusion coefficients to fix CMFD parameters to - dhat_reset : bool - Indicate whether :math:`\widehat{D}` nonlinear CMFD parameters should - be reset to zero before solving CMFD eigenproblem. display : dict Dictionary indicating which CMFD results to output. Note that CMFD k-effective will always be outputted. Acceptable keys are: @@ -312,7 +309,6 @@ class CMFDRun(object): self._tally_begin = 1 self._cmfd_begin = 1 self._ref_d = [] - self._dhat_reset = False self._display = {'balance': False, 'dominance': False, 'entropy': False, 'source': False} self._downscatter = False @@ -339,7 +335,6 @@ class CMFDRun(object): self._egrid = None self._albedo = None self._coremap = None - self._n_resets = 0 self._mesh_id = None self._tally_ids = None self._energy_filters = None @@ -425,10 +420,6 @@ class CMFDRun(object): def ref_d(self): return self._ref_d - @property - def dhat_reset(self): - return self._dhat_reset - @property def display(self): return self._display @@ -543,11 +534,6 @@ class CMFDRun(object): Iterable, Real) self._ref_d = np.array(diff_params) - @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('display', display, Mapping) @@ -762,22 +748,23 @@ class CMFDRun(object): calling :func:`openmc.capi.simulation_init` """ - # Configure CMFD parameters and tallies + # Configure CMFD parameters self._configure_cmfd() - # Initialize all arrays used for CMFD solver - self._allocate_cmfd() + # Create tally objects + self._create_cmfd_tally() - # Compute and store array indices used to build cross section - # arrays - self._precompute_array_indices() + if openmc.capi.master(): + # Compute and store array indices used to build cross section + # arrays + self._precompute_array_indices() - # Compute and store row and column indices used to build CMFD - # matrices - self._precompute_matrix_indices() + # Compute and store row and column indices used to build CMFD + # matrices + self._precompute_matrix_indices() - # Initialize all variables used for linear solver in C++ - self._initialize_linsolver() + # Initialize all variables used for linear solver in C++ + self._initialize_linsolver() # Initialize simulation openmc.capi.simulation_init() @@ -818,8 +805,9 @@ class CMFDRun(object): # Finalize simuation openmc.capi.simulation_finalize() - # Print out CMFD timing statistics - self._write_cmfd_timing_stats() + if openmc.capi.master(): + # Print out CMFD timing statistics + self._write_cmfd_timing_stats() def statepoint_write(self, filename=None): """Write all simulation parameters to statepoint @@ -943,53 +931,33 @@ class CMFDRun(object): def _write_cmfd_timing_stats(self): """Write CMFD timing stats to buffer after finalizing simulation""" - if openmc.capi.master(): - outstr = ("=====================> " - "CMFD TIMING STATISTICS <====================\n\n" - " Time in CMFD = {:.5E} seconds\n" - " Building matrices = {:.5E} seconds\n" - " Solving matrices = {:.5E} seconds\n") - print(outstr.format(self._time_cmfd, self._time_cmfdbuild, - self._time_cmfdsolve)) - sys.stdout.flush() + outstr = ("=====================> " + "CMFD TIMING STATISTICS <====================\n\n" + " Time in CMFD = {:.5E} seconds\n" + " Building matrices = {:.5E} seconds\n" + " Solving matrices = {:.5E} seconds\n") + print(outstr.format(self._time_cmfd, self._time_cmfdbuild, + self._time_cmfdsolve)) + sys.stdout.flush() def _configure_cmfd(self): """Initialize CMFD parameters and set CMFD input variables""" # Check if restarting simulation from statepoint file if not openmc.capi.settings.restart_run: - # Read in cmfd input defined in Python - self._read_cmfd_input() - - # Set up CMFD coremap - self._set_coremap() - - # Extract spatial and energy indices - nx, ny, nz, ng = self._indices - - # Allocate parameters that need to stored for tally window - self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0)) - self._flux_rate = np.zeros((nx, ny, nz, ng, 0)) - self._total_rate = np.zeros((nx, ny, nz, ng, 0)) - self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0)) - self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0)) - self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0)) - self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0)) - - # Initialize timers - self._time_cmfd = 0.0 - self._time_cmfdbuild = 0.0 - self._time_cmfdsolve = 0.0 - - # Initialize parameters for CMFD tally windows - self._set_tally_window() + # Define all variables necessary for running CMFD + self._initialize_cmfd() else: # Reset CMFD parameters from statepoint file path_statepoint = openmc.capi.settings.path_statepoint self._reset_cmfd(path_statepoint) - def _read_cmfd_input(self): - """Sets values of additional instance variables based on user input""" + def _initialize_cmfd(self): + """Sets values of CMFD instance variables based on user input, + separating between variables that only exist on all processes + and those that only exist on the master process + + """ # Print message to user and flush output to stdout if openmc.capi.settings.verbosity >= 7 and openmc.capi.master(): print(' Configuring CMFD parameters for simulation') @@ -1019,34 +987,55 @@ class CMFDRun(object): self._indices[3] = 1 self._energy_filters = False - # Set global albedo - if self._mesh.albedo is not None: - self._albedo = np.array(self._mesh.albedo) - else: - self._albedo = np.array([1., 1., 1., 1., 1., 1.]) - # Get acceleration map, otherwise set all regions to be accelerated if self._mesh.map is not None: check_length('CMFD coremap', self._mesh.map, np.product(self._indices[0:3])) - self._coremap = np.array(self._mesh.map) + if openmc.capi.master(): + self._coremap = np.array(self._mesh.map) else: - self._coremap = np.ones((np.product(self._indices[0:3])), - dtype=int) + if openmc.capi.master(): + self._coremap = np.ones((np.product(self._indices[0:3])), + dtype=int) # Check CMFD tallies accummulated before feedback turned on if self._feedback and self._cmfd_begin < self._tally_begin: raise ValueError('Tally begin must be less than or equal to ' 'CMFD begin') - # Set number of batches where cmfd tallies should be reset - self._n_resets = len(self._reset) + # Initialize parameters for CMFD tally windows + self._set_tally_window() - # Create tally objects - self._create_cmfd_tally() + # Define all variables that will exist only on master process + if openmc.capi.master(): + # Set global albedo + if self._mesh.albedo is not None: + self._albedo = np.array(self._mesh.albedo) + else: + self._albedo = np.array([1., 1., 1., 1., 1., 1.]) + + # Set up CMFD coremap + self._set_coremap() + + # Extract spatial and energy indices + nx, ny, nz, ng = self._indices + + # Allocate parameters that need to be stored for tally window + self._openmc_src_rate = np.zeros((nx, ny, nz, ng, 0)) + self._flux_rate = np.zeros((nx, ny, nz, ng, 0)) + self._total_rate = np.zeros((nx, ny, nz, ng, 0)) + self._p1scatt_rate = np.zeros((nx, ny, nz, ng, 0)) + self._scatt_rate = np.zeros((nx, ny, nz, ng, ng, 0)) + self._nfiss_rate = np.zeros((nx, ny, nz, ng, ng, 0)) + self._current_rate = np.zeros((nx, ny, nz, 12, ng, 0)) + + # Initialize timers + self._time_cmfd = 0.0 + self._time_cmfdbuild = 0.0 + self._time_cmfdsolve = 0.0 def _reset_cmfd(self, filename): - """Reset all CMFD parameters from statepoint + """Reset all CMFD parameters from statepoint Parameters ---------- @@ -1065,32 +1054,28 @@ class CMFDRun(object): print(' Loading CMFD data from {}...'.format(filename)) sys.stdout.flush() cmfd_group = f['cmfd'] + + # Define variables that exist on all processes self._cmfd_on = cmfd_group.attrs['cmfd_on'] self._feedback = cmfd_group.attrs['feedback'] self._cmfd_begin = cmfd_group.attrs['cmfd_begin'] self._tally_begin = cmfd_group.attrs['tally_begin'] - self._time_cmfd = cmfd_group.attrs['time_cmfd'] - self._time_cmfdbuild = cmfd_group.attrs['time_cmfdbuild'] - self._time_cmfdsolve = cmfd_group.attrs['time_cmfdsolve'] - self._window_size = cmfd_group.attrs['window_size'] - self._window_type = cmfd_group.attrs['window_type'] self._k_cmfd = list(cmfd_group['k_cmfd']) self._dom = list(cmfd_group['dom']) self._src_cmp = list(cmfd_group['src_cmp']) self._balance = list(cmfd_group['balance']) self._entropy = list(cmfd_group['entropy']) self._reset = list(cmfd_group['reset']) - self._albedo = cmfd_group['albedo'][()] - self._coremap = cmfd_group['coremap'][()] self._egrid = cmfd_group['egrid'][()] self._indices = cmfd_group['indices'][()] - self._current_rate = cmfd_group['current_rate'][()] - self._flux_rate = cmfd_group['flux_rate'][()] - self._nfiss_rate = cmfd_group['nfiss_rate'][()] - self._openmc_src_rate = cmfd_group['openmc_src_rate'][()] - self._p1scatt_rate = cmfd_group['p1scatt_rate'][()] - self._scatt_rate = cmfd_group['scatt_rate'][()] - self._total_rate = cmfd_group['total_rate'][()] + default_egrid = np.array([_ENERGY_MIN_NEUTRON, + _ENERGY_MAX_NEUTRON]) + self._energy_filters = not np.array_equal(self._egrid, + default_egrid) + self._window_size = cmfd_group.attrs['window_size'] + self._window_type = cmfd_group.attrs['window_type'] + self._reset_every = (self._window_type == 'expanding' or + self._window_type == 'rolling') # Overwrite CMFD mesh properties cmfd_mesh_name = 'mesh ' + str(cmfd_group.attrs['mesh_id']) @@ -1100,49 +1085,21 @@ class CMFDRun(object): self._mesh.upper_right = cmfd_mesh['upper_right'][()] self._mesh.width = cmfd_mesh['width'][()] - # Store tally ids from statepoint run - sp_tally_ids = list(cmfd_group['tally_ids']) - - # Set CMFD variables not in statepoint file - default_egrid = np.array([_ENERGY_MIN_NEUTRON, _ENERGY_MAX_NEUTRON]) - self._energy_filters = not np.array_equal(self._egrid, default_egrid) - self._n_resets = len(self._reset) - self._mat_dim = np.max(self._coremap) + 1 - self._reset_every = (self._window_type == 'expanding' or - self._window_type == 'rolling') - - # Recreate CMFD tallies in memory - self._create_cmfd_tally() - - def _allocate_cmfd(self): - """Allocates all numpy arrays and lists used in CMFD algorithm""" - # Extract spatial and energy indices - nx, ny, nz, ng = self._indices - - # Allocate dimensions for each mesh cell - self._hxyz = np.zeros((nx, ny, nz, 3)) - self._hxyz[:] = openmc.capi.meshes[self._mesh_id].width - - # Allocate flux, cross sections and diffusion coefficient - self._flux = np.zeros((nx, ny, nz, ng)) - self._totalxs = np.zeros((nx, ny, nz, ng)) - self._p1scattxs = np.zeros((nx, ny, nz, ng)) - self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing - self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing - self._diffcof = np.zeros((nx, ny, nz, ng)) - - # Allocate dtilde and dhat - self._dtilde = np.zeros((nx, ny, nz, ng, 6)) - self._dhat = np.zeros((nx, ny, nz, ng, 6)) - - # Set reference diffusion parameters - if list(self._ref_d): - 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]: - raise OpenMCError('Number of reference diffusion parameters ' - 'must equal number of CMFD energy groups') + # Define variables that exist only on master process + if openmc.capi.master(): + self._time_cmfd = cmfd_group.attrs['time_cmfd'] + self._time_cmfdbuild = cmfd_group.attrs['time_cmfdbuild'] + self._time_cmfdsolve = cmfd_group.attrs['time_cmfdsolve'] + self._albedo = cmfd_group['albedo'][()] + self._coremap = cmfd_group['coremap'][()] + self._current_rate = cmfd_group['current_rate'][()] + self._flux_rate = cmfd_group['flux_rate'][()] + self._nfiss_rate = cmfd_group['nfiss_rate'][()] + self._openmc_src_rate = cmfd_group['openmc_src_rate'][()] + self._p1scatt_rate = cmfd_group['p1scatt_rate'][()] + self._scatt_rate = cmfd_group['scatt_rate'][()] + self._total_rate = cmfd_group['total_rate'][()] + self._mat_dim = np.max(self._coremap) + 1 def _set_tally_window(self): """Sets parameters to handle different tally window options""" @@ -1169,7 +1126,7 @@ class CMFDRun(object): self._cmfd_on = True # Check to reset tallies - if ((self._n_resets > 0 and current_batch in self._reset) + if ((len(self._reset) > 0 and current_batch in self._reset) or self._reset_every): self._cmfd_tally_reset() @@ -2206,12 +2163,37 @@ class CMFDRun(object): (ng * num_accel))) def _precompute_array_indices(self): - """Computes the indices used to populate certain cross section arrays. - These indices are used in _compute_dtilde and _compute_dhat + """Initializes cross section arrays and computes the indices + used to populate dtilde and dhat """ # Extract spatial indices - nx, ny, nz = self._indices[:3] + nx, ny, nz, ng = self._indices + + # Allocate dimensions for each mesh cell + self._hxyz = np.zeros((nx, ny, nz, 3)) + self._hxyz[:] = openmc.capi.meshes[self._mesh_id].width + + # Allocate flux, cross sections and diffusion coefficient + self._flux = np.zeros((nx, ny, nz, ng)) + self._totalxs = np.zeros((nx, ny, nz, ng)) + self._p1scattxs = np.zeros((nx, ny, nz, ng)) + self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing + self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Incoming, outgoing + self._diffcof = np.zeros((nx, ny, nz, ng)) + + # Allocate dtilde and dhat + self._dtilde = np.zeros((nx, ny, nz, ng, 6)) + self._dhat = np.zeros((nx, ny, nz, ng, 6)) + + # Set reference diffusion parameters + if list(self._ref_d): + 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]: + raise OpenMCError('Number of reference diffusion parameters ' + 'must equal number of CMFD energy groups') # Logical for determining whether region of interest is accelerated # region From 603ea4ec14a6d71d0dbfa7afe061277a8fadc473 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 1 Aug 2019 12:27:33 -0400 Subject: [PATCH 12/16] Free CMFD memory on master process --- src/finalize.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/finalize.cpp b/src/finalize.cpp index 633cf2991d..4d93c6c66c 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -44,7 +44,9 @@ void free_memory() free_memory_mesh(); free_memory_tally(); free_memory_bank(); - free_memory_cmfd(); + if (mpi::master) { + free_memory_cmfd(); + } #ifdef DAGMC free_memory_dagmc(); #endif From 0ce75aa21b54bf350fb1869a61501b39eb21a8d0 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Wed, 25 Sep 2019 18:29:48 -0400 Subject: [PATCH 13/16] Address Paul comments except setting threads as user param --- include/openmc/capi.h | 2 +- openmc/capi/core.py | 2 +- openmc/cmfd.py | 73 +++++++++++-------- src/cmfd_solver.cpp | 32 ++++++-- tests/regression_tests/cmfd_feed/test.py | 8 +- tests/regression_tests/cmfd_feed_2g/test.py | 2 +- .../cmfd_feed_expanding_window/test.py | 2 +- tests/regression_tests/cmfd_feed_ng/test.py | 2 +- .../regression_tests/cmfd_feed_ref_d/test.py | 2 +- .../cmfd_feed_rolling_window/test.py | 2 +- tests/regression_tests/cmfd_nofeed/test.py | 2 +- tests/regression_tests/cmfd_restart/test.py | 4 +- 12 files changed, 84 insertions(+), 49 deletions(-) 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] From 8e0b6160751398985d464e9e7b843b4d803895ed Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 4 Nov 2019 22:35:05 -0500 Subject: [PATCH 14/16] Create use_all_threads variable instead of n_threads --- include/openmc/capi.h | 2 +- openmc/capi/core.py | 2 +- openmc/cmfd.py | 23 +++++++++++------------ src/cmfd_solver.cpp | 18 ++++++------------ 4 files changed, 19 insertions(+), 26 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 0acc29e780..dbdf713973 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, int n_threads); + const int* map, bool use_all_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 fe911d3a52..546121c98c 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_int] + c_double, _array_1d_int, _array_1d_int, c_bool] _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 c711c9b468..445eee7f4a 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -295,6 +295,8 @@ class CMFDRun(object): Time for building CMFD matrices, in seconds time_cmfdsolve : float Time for solving CMFD matrix equations, in seconds + use_all_threads : bool + Whether to use all threads allocated to OpenMC for CMFD solver intracomm : mpi4py.MPI.Intracomm or None MPI intercommunicator for running MPI commands @@ -328,7 +330,7 @@ class CMFDRun(object): self._window_type = 'none' self._window_size = 10 self._intracomm = None - self._n_threads = 1 + self._use_all_threads = False # External variables used during runtime but users cannot control self._set_reference_params = False @@ -494,8 +496,8 @@ class CMFDRun(object): return self._indices @property - def n_threads(self): - return self._n_threads + def use_all_threads(self): + return self._use_all_threads @property def cmfd_src(self): @@ -683,11 +685,10 @@ 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 + @use_all_threads.setter + def use_all_threads(self, use_all_threads): + check_type('CMFD use all threads', use_all_threads, bool) + self._use_all_threads = use_all_threads def run(self, **kwargs): """Run OpenMC with coarse mesh finite difference acceleration @@ -704,8 +705,7 @@ class CMFDRun(object): """ with self.run_in_memory(**kwargs): for _ in self.iter_batches(): - print('done') - #pass + pass @contextmanager def run_in_memory(self, **kwargs): @@ -800,7 +800,6 @@ class CMFDRun(object): # Run next batch status = openmc.capi.next_batch() - print('2') # Perform CMFD calculations self._execute_cmfd() @@ -918,7 +917,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._n_threads + self._spectral, self._indices, coremap, self._use_all_threads return openmc.capi._dll.openmc_initialize_linsolver(*args) def _write_cmfd_output(self): diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 40be440d5f..6c432f595c 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -32,9 +32,7 @@ int nx, ny, nz, ng; xt::xtensor indexmap; -int n_threads; - -int n_threads_reset; +int use_all_threads; } // namespace cmfd @@ -108,7 +106,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) num_threads(cmfd::n_threads) + #pragma omp parallel for reduction (+:err) if(cmfd::use_all_threads) for (int irow = 0; irow < cmfd::dim; irow++) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -175,7 +173,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) num_threads(cmfd::n_threads) + #pragma omp parallel for reduction (+:err) if(cmfd::use_all_threads) for (int irow = 0; irow < cmfd::dim; irow+=2) { int g, i, j, k; matrix_to_indices(irow, g, i, j, k); @@ -311,7 +309,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, int n_threads) + const int* map, bool use_all_threads) { // Store elements of indptr for (int i = 0; i < len_indptr; i++) @@ -339,12 +337,8 @@ void openmc_initialize_linsolver(const int* indptr, int len_indptr, 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 + // Use all threads allocated to OpenMC simulation to run CMFD solver + cmfd::use_all_threads = use_all_threads; } //============================================================================== From 62c50ea302ecc111f1c53fe1891c85fe92628a86 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 4 Nov 2019 23:05:34 -0500 Subject: [PATCH 15/16] Add test for multithreaded CMFD --- src/cmfd_solver.cpp | 12 +++--------- tests/regression_tests/cmfd_feed/test.py | 24 ++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 6c432f595c..1fdbdf5668 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -350,20 +350,14 @@ extern "C" int openmc_run_linsolver(const double* A_data, const double* b, double* x, double tol) { - int result; - switch (cmfd::ng) { case 1: - result = cmfd_linsolver_1g(A_data, b, x, tol); - break; + return cmfd_linsolver_1g(A_data, b, x, tol); case 2: - result = cmfd_linsolver_2g(A_data, b, x, tol); - break; + return cmfd_linsolver_2g(A_data, b, x, tol); default: - result = cmfd_linsolver_ng(A_data, b, x, tol); - break; + return cmfd_linsolver_ng(A_data, b, x, tol); } - 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 6b56635f4a..b4dbe682b8 100644 --- a/tests/regression_tests/cmfd_feed/test.py +++ b/tests/regression_tests/cmfd_feed/test.py @@ -140,3 +140,27 @@ def test_cmfd_feed(): # Initialize and run CMFD test harness harness = CMFDTestHarness('statepoint.20.h5', cmfd_run) harness.main() + +def test_cmfd_multithread(): + """Test 1 group CMFD solver with CMFD feedback""" + # Initialize and set CMFD mesh + cmfd_mesh = cmfd.CMFDMesh() + cmfd_mesh.lower_left = (-10.0, -1.0, -1.0) + cmfd_mesh.upper_right = (10.0, 1.0, 1.0) + cmfd_mesh.dimension = (10, 1, 1) + cmfd_mesh.albedo = (0.0, 0.0, 1.0, 1.0, 1.0, 1.0) + + # Initialize and run CMFDRun object + cmfd_run = cmfd.CMFDRun() + cmfd_run.mesh = cmfd_mesh + cmfd_run.tally_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] + cmfd_run.use_all_threads = True + cmfd_run.run() + + # Initialize and run CMFD test harness + harness = CMFDTestHarness('statepoint.20.h5', cmfd_run) + harness.main() From 6367052fa55298e4287478703e5d61340e735ffa Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Mon, 4 Nov 2019 23:07:14 -0500 Subject: [PATCH 16/16] Update description for regression test --- tests/regression_tests/cmfd_feed/test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/cmfd_feed/test.py b/tests/regression_tests/cmfd_feed/test.py index b4dbe682b8..7723a47b1f 100644 --- a/tests/regression_tests/cmfd_feed/test.py +++ b/tests/regression_tests/cmfd_feed/test.py @@ -142,7 +142,7 @@ def test_cmfd_feed(): harness.main() def test_cmfd_multithread(): - """Test 1 group CMFD solver with CMFD feedback""" + """Test 1 group CMFD solver with all available threads""" # Initialize and set CMFD mesh cmfd_mesh = cmfd.CMFDMesh() cmfd_mesh.lower_left = (-10.0, -1.0, -1.0)