From 0da3bf338180255d80a3a997bebe8445d4667a82 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 9 Jul 2019 19:54:40 -0400 Subject: [PATCH 01/24] 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 f7ef57546..c46ace3f7 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/24] 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 a37a78fa9..439dcde6a 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 c46ace3f7..472107d0c 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 07ade9c70..e831d8215 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/24] 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 742cdeabf..ca52502c3 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/24] 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 472107d0c..11d6d97b5 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/24] 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 11d6d97b5..27d594390 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/24] 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 ca52502c3..8fa406bc5 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 906c631ec..87370a894 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 ba4d21609..8f0fcda68 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 b67bef0f7..98b7a5e82 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 e8671e95c..35cf3c289 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 dce674f90..6b9ff60c5 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 06440d81f..4a26ded63 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 d033ab8dd..925fec32d 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 518052683..98a8f8d06 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 802b5b700..07eb0741b 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 7d0895c2f..35c18a38f 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 7ed92d024..369ef0349 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/24] 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 8fa406bc5..ac425d283 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/24] 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 ac425d283..0a18dffdf 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/24] 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 0a18dffdf..e4976c92c 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/24] 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 e4976c92c..b0a4f55b1 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/24] 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 b0a4f55b1..69a24ba61 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/24] 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 633cf2991..4d93c6c66 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/24] 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 3d2e8f57b..0acc29e78 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 a470f0665..fe911d3a5 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 69a24ba61..c711c9b46 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 27d594390..40be440d5 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 87370a894..6b56635f4 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 8f0fcda68..d3af8998b 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 35cf3c289..964d4f225 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 6b9ff60c5..a2a522e9c 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 925fec32d..120d94b6b 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 07eb0741b..2c7b7f242 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 35c18a38f..7ab72f86a 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 369ef0349..8f8410a7a 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/24] 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 0acc29e78..dbdf71397 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 fe911d3a5..546121c98 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 c711c9b46..445eee7f4 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 40be440d5..6c432f595 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/24] 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 6c432f595..1fdbdf566 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 6b56635f4..b4dbe682b 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/24] 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 b4dbe682b..7723a47b1 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) From df5727cf3e0e3e63603193df147342bf9ea922e9 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sat, 14 Dec 2019 21:02:48 -0500 Subject: [PATCH 17/24] Fix setting multiplicity matrix for histogram MGMC XS --- openmc/mgxs_library.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index b98902030..34c36ec17 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1638,9 +1638,9 @@ class XSdata(object): xs_type=xs_type, moment=0, subdomains=subdomain) if scatter.scatter_format == 'histogram': - scatt = np.sum(scatt, axis=0) + scatt = np.sum(scatt, axis=2) if nuscatter.scatter_format == 'histogram': - nuscatt = np.sum(nuscatt, axis=0) + nuscatt = np.sum(nuscatt, axis=2) self._multiplicity_matrix[i] = np.divide(nuscatt, scatt) self._multiplicity_matrix[i] = \ From e97a69069ebbf76cec325620843a49929dbaaa6d Mon Sep 17 00:00:00 2001 From: guillaume Date: Sat, 14 Dec 2019 21:17:00 -0500 Subject: [PATCH 18/24] Handle case where user forgot to tally multiplicities or nu-scatter more gracefully --- openmc/mgxs/library.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 83d3178ff..352c5b173 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -1167,6 +1167,16 @@ class Library(object): np.subtract(xsdata._total[i], np.sum(np.sum( xsdata._scatter_matrix[i][:, :, :, :, :], axis=4), axis=3)) + # if only scatter matrices have been tallied, multiplicity cannot + # be accounted for + else: + msg = 'Scatter multiplicity (such as (n,xn) reactions) ' \ + 'are ignored since multiplicity or nu-scatter matrices '\ + 'were not tallied for '+xsdata_name + warn(msg, RuntimeWarning) + xsdata.set_scatter_matrix_mgxs(scatt_mgxs, xs_type=xs_type, + nuclide=[nuclide], + subdomain=subdomain) return xsdata From 80dc9967179e5152e4b9a1b5758baf51d147e5ec Mon Sep 17 00:00:00 2001 From: guillaume Date: Sat, 14 Dec 2019 21:27:41 -0500 Subject: [PATCH 19/24] Check for multiplicity matrix tallies in routine meant to check libraries for MGMC library creation --- openmc/mgxs/library.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 352c5b173..0f10826eb 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -1420,6 +1420,8 @@ class Library(object): fixed source problem could be the target. - Fission and kappa-fission are not required as they are only needed to support tallies the user may wish to request. + - Scattering multiplicity should have been tallied, either using a + multiplicity or scatter and nu-scatter matrix tally. See also -------- @@ -1470,6 +1472,16 @@ class Library(object): '"scatter matrix", or "consistent scatter matrix" MGXS ' 'type is required.') + # Make sure there is some kind of a scattering multiplicity matrix data + if 'multiplicity matrix' not in self.mgxs_types and \ + ('scatter matrix' not in self.mgxs_types or + 'nu-scatter matrix' not in self.mgxs_types) and\ + ('consistent scatter matrix' not in self.mgxs_types or + 'consistent nu-scatter matrix' not in self.mgxs_types): + error_flag = True + warn('A "multiplicity matrix" or both a "scatter" and "nu-scatter" ' + 'matrix MGXS type(s) is/are required.') + # Ensure absorption is present if 'absorption' not in self.mgxs_types: error_flag = True From d9195b93f8f2bf5876370886d3c3678a59015c90 Mon Sep 17 00:00:00 2001 From: guillaume Date: Mon, 16 Dec 2019 17:35:19 -0500 Subject: [PATCH 20/24] Address @nelsonag 's review --- openmc/mgxs/library.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 0f10826eb..9b62634c7 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -1172,7 +1172,7 @@ class Library(object): else: msg = 'Scatter multiplicity (such as (n,xn) reactions) ' \ 'are ignored since multiplicity or nu-scatter matrices '\ - 'were not tallied for '+xsdata_name + 'were not tallied for ' + xsdata_name warn(msg, RuntimeWarning) xsdata.set_scatter_matrix_mgxs(scatt_mgxs, xs_type=xs_type, nuclide=[nuclide], @@ -1420,8 +1420,9 @@ class Library(object): fixed source problem could be the target. - Fission and kappa-fission are not required as they are only needed to support tallies the user may wish to request. - - Scattering multiplicity should have been tallied, either using a - multiplicity or scatter and nu-scatter matrix tally. + - Scattering multiplicity should have been tallied for increased model + accuracy, either using a multiplicity or scatter and nu-scatter matrix + tally. See also -------- @@ -1478,9 +1479,8 @@ class Library(object): 'nu-scatter matrix' not in self.mgxs_types) and\ ('consistent scatter matrix' not in self.mgxs_types or 'consistent nu-scatter matrix' not in self.mgxs_types): - error_flag = True warn('A "multiplicity matrix" or both a "scatter" and "nu-scatter" ' - 'matrix MGXS type(s) is/are required.') + 'matrix MGXS type(s) should be provided.') # Ensure absorption is present if 'absorption' not in self.mgxs_types: From a50a342f273b0f4b4818c37ec8a5f83af8534209 Mon Sep 17 00:00:00 2001 From: guillaume Date: Mon, 16 Dec 2019 17:38:29 -0500 Subject: [PATCH 21/24] Add MGMC loading of regular transport cross sections --- openmc/mgxs/library.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 9b62634c7..701f7d8eb 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -999,6 +999,11 @@ class Library(object): xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], subdomain=subdomain) + elif 'transport' in self.mgxs_types and self.correction == 'P0': + mymgxs = self.get_mgxs(domain, 'transport') + xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], + subdomain=subdomain) + elif 'total' in self.mgxs_types: mymgxs = self.get_mgxs(domain, 'total') xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide], @@ -1170,7 +1175,7 @@ class Library(object): # if only scatter matrices have been tallied, multiplicity cannot # be accounted for else: - msg = 'Scatter multiplicity (such as (n,xn) reactions) ' \ + msg = 'Scatter multiplicity (such as (n,xn) reactions) '\ 'are ignored since multiplicity or nu-scatter matrices '\ 'were not tallied for ' + xsdata_name warn(msg, RuntimeWarning) From aa768abbc3c730c58613dee72110f4e38162dcc5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jan 2020 16:33:14 -0600 Subject: [PATCH 22/24] Happy new year 2020! --- LICENSE | 2 +- docs/source/conf.py | 2 +- docs/source/license.rst | 2 +- man/man1/openmc.1 | 2 +- src/output.cpp | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/LICENSE b/LICENSE index 328f3ee46..c145d7a0d 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011-2019 Massachusetts Institute of Technology and OpenMC contributors +Copyright (c) 2011-2020 Massachusetts Institute of Technology and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/docs/source/conf.py b/docs/source/conf.py index 68caff888..ebda1bdad 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -72,7 +72,7 @@ master_doc = 'index' # General information about the project. project = 'OpenMC' -copyright = '2011-2019, Massachusetts Institute of Technology and OpenMC contributors' +copyright = '2011-2020, Massachusetts Institute of Technology and OpenMC contributors' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the diff --git a/docs/source/license.rst b/docs/source/license.rst index 98f6655bc..e9b71d882 100644 --- a/docs/source/license.rst +++ b/docs/source/license.rst @@ -4,7 +4,7 @@ License Agreement ================= -Copyright © 2011-2019 Massachusetts Institute of Technology and OpenMC contributors +Copyright © 2011-2020 Massachusetts Institute of Technology and OpenMC contributors Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in diff --git a/man/man1/openmc.1 b/man/man1/openmc.1 index 229b64fd1..168fab686 100644 --- a/man/man1/openmc.1 +++ b/man/man1/openmc.1 @@ -54,7 +54,7 @@ Indicates the default path to an HDF5 file that contains multi-group cross section libraries if the user has not specified the tag in .I materials.xml\fP. .SH LICENSE -Copyright \(co 2011-2019 Massachusetts Institute of Technology and OpenMC +Copyright \(co 2011-2020 Massachusetts Institute of Technology and OpenMC contributors. .PP Permission is hereby granted, free of charge, to any person obtaining a copy of diff --git a/src/output.cpp b/src/output.cpp index a388b15d7..c8b91c4dc 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -72,7 +72,7 @@ void title() // Write version information std::cout << " | The OpenMC Monte Carlo Code\n" << - " Copyright | 2011-2019 MIT and OpenMC contributors\n" << + " Copyright | 2011-2020 MIT and OpenMC contributors\n" << " License | http://openmc.readthedocs.io/en/latest/license.html\n" << " Version | " << VERSION_MAJOR << '.' << VERSION_MINOR << '.' << VERSION_RELEASE << (VERSION_DEV ? "-dev" : "") << '\n'; From 2b1380519d8c7614fd74aa56cd66f8dac6637c04 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 20 Dec 2019 15:55:21 -0500 Subject: [PATCH 23/24] Adding data_type kwarg to DataLibrary.get_by_material Co-Authored-By: Paul Romano Update tests/unit_tests/test_data_misc.py Co-Authored-By: Paul Romano --- openmc/data/library.py | 7 +++++-- openmc/plotter.py | 8 +++++--- tests/unit_tests/test_data_misc.py | 2 +- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index 2339ca1c9..3830334d2 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -23,13 +23,16 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] - def get_by_material(self, name): + def get_by_material(self, name, data_type='neutron'): """Return the library dictionary containing a given material. Parameters ---------- name : str Name of material, e.g. 'Am241' + data_type : str + Name of data type, e.g. 'neutron', 'photon', 'wmp', + or 'thermal' Returns ------- @@ -39,7 +42,7 @@ class DataLibrary(EqualityMixin): """ for library in self.libraries: - if name in library['materials']: + if name in library['materials'] and data_type in library['type']: return library return None diff --git a/openmc/plotter.py b/openmc/plotter.py index af317dfdf..33f1e4797 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -562,14 +562,16 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., if isinstance(this, openmc.Material): for sab_name in this._sab: sab = openmc.data.ThermalScattering.from_hdf5( - library.get_by_material(sab_name)['path']) + library.get_by_material(sab_name, data_type='thermal')['path']) for nuc in sab.nuclides: - sabs[nuc] = library.get_by_material(sab_name)['path'] + sabs[nuc] = library.get_by_material(sab_name, + data_type='thermal')['path'] else: if sab_name: sab = openmc.data.ThermalScattering.from_hdf5(sab_name) for nuc in sab.nuclides: - sabs[nuc] = library.get_by_material(sab_name)['path'] + sabs[nuc] = library.get_by_material(sab_name, + data_type='thermal')['path'] # Now we can create the data sets to be plotted xs = {} diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index 2585e9315..8244ca7f0 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -18,7 +18,7 @@ def test_data_library(tmpdir): assert f['type'] == 'neutron' assert 'U235' in f['materials'] - f = lib.get_by_material('c_H_in_H2O') + f = lib.get_by_material('c_H_in_H2O', data_type='thermal') assert f['type'] == 'thermal' assert 'c_H_in_H2O' in f['materials'] From 96c7427eafb951dbe1cb8fd539cc403e4db7a67b Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 31 Dec 2019 12:11:49 -0500 Subject: [PATCH 24/24] Adding CylindricalIndependent source type ammended docstrings to mention local reference frame --- docs/source/io_formats/settings.rst | 53 ++++--- docs/source/pythonapi/stats.rst | 1 + include/openmc/distribution_spatial.h | 20 +++ openmc/stats/multivariate.py | 149 ++++++++++++++++-- src/distribution_spatial.cpp | 66 ++++++++ src/source.cpp | 2 + tests/regression_tests/source/inputs_true.dat | 13 ++ .../regression_tests/source/results_true.dat | 2 +- tests/regression_tests/source/test.py | 8 +- 9 files changed, 277 insertions(+), 37 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index c14767e46..3a5e347f9 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -423,16 +423,19 @@ attributes/sub-elements: :type: The type of spatial distribution. Valid options are "box", "fission", - "point", "cartesian", and "spherical". A "box" spatial distribution has - coordinates sampled uniformly in a parallelepiped. A "fission" spatial - distribution samples locations from a "box" distribution but only - locations in fissionable materials are accepted. A "point" spatial - distribution has coordinates specified by a triplet. An "cartesian" - spatial distribution specifies independent distributions of x-, y-, and - z-coordinates. A "spherical" spatial distribution specifies independent - distributions of r-, theta-, and phi-coordinates where theta is the angle - with respect to the z-axis, phi is the azimuthal angle, and the sphere is - centered on the coordinate (x0,y0,z0). + "point", "cartesian", "cylindrical", and "spherical". A "box" spatial + distribution has coordinates sampled uniformly in a parallelepiped. A + "fission" spatial distribution samples locations from a "box" + distribution but only locations in fissionable materials are accepted. + A "point" spatial distribution has coordinates specified by a triplet. + A "cartesian" spatial distribution specifies independent distributions of + x-, y-, and z-coordinates. A "cylindrical" spatial distribution specifies + independent distributions of r-, phi-, and z-coordinates where phi is the + azimuthal angle and the origin for the cylindrical coordinate system is + specified by origin. A "spherical" spatial distribution specifies + independent distributions of r-, theta-, and phi-coordinates where theta + is the angle with respect to the z-axis, phi is the azimuthal angle, and + the sphere is centered on the coordinate (x0,y0,z0). *Default*: None @@ -449,6 +452,9 @@ attributes/sub-elements: For an "cartesian" distribution, no parameters are specified. Instead, the ``x``, ``y``, and ``z`` elements must be specified. + + For a "cylindrical" distribution, no parameters are specified. Instead, + the ``r``, ``phi``, ``z``, and ``origin`` elements must be specified. For a "spherical" distribution, no parameters are specified. Instead, the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified. @@ -468,15 +474,16 @@ attributes/sub-elements: :ref:`univariate`). :z: - For an "cartesian" distribution, this element specifies the distribution - of z-coordinates. The necessary sub-elements/attributes are those of a - univariate probability distribution (see the description in - :ref:`univariate`). + For both "cartesian" and "cylindrical" distributions, this element + specifies the distribution of z-coordinates. The necessary + sub-elements/attributes are those of a univariate probability + distribution (see the description in :ref:`univariate`). :r: - For a "spherical" distribution, this element specifies the distribution - of r-coordinates. The necessary sub-elements/attributes are those of a - univariate probability distribution (see the description in + For "cylindrical" and "spherical" distributions, this element specifies + the distribution of r-coordinates (cylindrical radius and spherical + radius, respectively). The necessary sub-elements/attributes are those + of a univariate probability distribution (see the description in :ref:`univariate`). :theta: @@ -486,14 +493,14 @@ attributes/sub-elements: :ref:`univariate`). :phi: - For a "spherical" distribution, this element specifies the distribution - of phi-coordinates. The necessary sub-elements/attributes are those of a - univariate probability distribution (see the description in - :ref:`univariate`). + For "cylindrical" and "spherical" distributions, this element specifies + the distribution of phi-coordinates. The necessary + sub-elements/attributes are those of a univariate probability + distribution (see the description in :ref:`univariate`). :origin: - For a "spherical" distribution, this element specifies the coordinates of - the center of the sphere. + For "cylindrical and "spherical" distributions, this element specifies + the coordinates for the origin of the coordinate system. :angle: An element specifying the angular distribution of source sites. This element diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index cacbab31b..d14aff558 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -46,6 +46,7 @@ Spatial Distributions openmc.stats.Spatial openmc.stats.CartesianIndependent + openmc.stats.CylindricalIndependent openmc.stats.SphericalIndependent openmc.stats.Box openmc.stats.Point diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 53857fac3..f3c7bb67d 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -38,6 +38,26 @@ private: UPtrDist z_; //!< Distribution of z coordinates }; +//============================================================================== +//! Distribution of points specified by cylindrical coordinates r,phi,z +//============================================================================== + +class CylindricalIndependent : public SpatialDistribution { +public: + explicit CylindricalIndependent(pugi::xml_node node); + + //! Sample a position from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled position + Position sample(uint64_t* seed) const; +private: + UPtrDist r_; //!< Distribution of r coordinates + UPtrDist phi_; //!< Distribution of phi coordinates + UPtrDist z_; //!< Distribution of z coordinates + Position origin_; //!< Cartesian coordinates of the cylinder center +}; + + //============================================================================== //! Distribution of points specified by spherical coordinates r,theta,phi //============================================================================== diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 6a97e6192..3d258ee35 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -272,6 +272,8 @@ class Spatial(metaclass=ABCMeta): distribution = get_text(elem, 'type') if distribution == 'cartesian': return CartesianIndependent.from_xml_element(elem) + elif distribution == 'cylindrical': + return CylindricalIndependent.from_xml_element(elem) elif distribution == 'spherical': return SphericalIndependent.from_xml_element(elem) elif distribution == 'box' or distribution == 'fission': @@ -377,7 +379,7 @@ class CartesianIndependent(Spatial): class SphericalIndependent(Spatial): - """Spatial distribution represented in spherical coordinates. + r"""Spatial distribution represented in spherical coordinates. This distribution allows one to specify coordinates whose :math:`r`, :math:`\theta`, and :math:`\phi` components are sampled independently from @@ -386,26 +388,31 @@ class SphericalIndependent(Spatial): Parameters ---------- r : openmc.stats.Univariate - Distribution of r-coordinates + Distribution of r-coordinates in a reference frame specified by + the origin parameter theta : openmc.stats.Univariate - Distribution of theta-coordinates (angle relative to the z-axis) + Distribution of theta-coordinates (angle relative to the z-axis) in a + reference frame specified by the origin parameter phi : openmc.stats.Univariate - Distribution of phi-coordinates (azimuthal angle) + Distribution of phi-coordinates (azimuthal angle) in a reference frame + specified by the origin parameter origin: Iterable of float, optional - coordinates (x0, y0, z0) of the center of the sphere. Defaults to - (0.0, 0.0, 0.0) + coordinates (x0, y0, z0) of the center of the spherical reference frame + for the source. Defaults to (0.0, 0.0, 0.0) Attributes ---------- r : openmc.stats.Univariate - Distribution of r-coordinates + Distribution of r-coordinates in the local reference frame theta : openmc.stats.Univariate - Distribution of theta-coordinates (angle relative to the z-axis) + Distribution of theta-coordinates (angle relative to the z-axis) in the + local reference frame phi : openmc.stats.Univariate - Distribution of phi-coordinates (azimuthal angle) + Distribution of phi-coordinates (azimuthal angle) in the local + reference frame origin: Iterable of float, optional - coordinates (x0, y0, z0) of the center of the sphere. Defaults to - (0.0, 0.0, 0.0) + coordinates (x0, y0, z0) of the center of the spherical reference + frame. Defaults to (0.0, 0.0, 0.0) """ @@ -491,6 +498,126 @@ class SphericalIndependent(Spatial): origin = [float(x) for x in elem.get('origin').split()] return cls(r, theta, phi, origin=origin) +class CylindricalIndependent(Spatial): + """Spatial distribution represented in cylindrical coordinates. + + This distribution allows one to specify coordinates whose :math:`r`, + :math:`\phi`, and :math:`z` components are sampled independently from + one another and in a reference frame whose origin is specified by the + coordinates (x0, y0, z0). + + Parameters + ---------- + r : openmc.stats.Univariate + Distribution of r-coordinates in a reference frame specified by the + origin parameter + phi : openmc.stats.Univariate + Distribution of phi-coordinates (azimuthal angle) in a reference frame + specified by the origin parameter + z : openmc.stats.Univariate + Distribution of z-coordinates in a reference frame specified by the + origin parameter + origin: Iterable of float, optional + coordinates (x0, y0, z0) of the center of the cylindrical reference + frame. Defaults to (0.0, 0.0, 0.0) + + Attributes + ---------- + r : openmc.stats.Univariate + Distribution of r-coordinates in the local reference frame + phi : openmc.stats.Univariate + Distribution of phi-coordinates (azimuthal angle) in the local + reference frame + z : openmc.stats.Univariate + Distribution of z-coordinates in the local reference frame + origin: Iterable of float, optional + coordinates (x0, y0, z0) of the center of the cylindrical reference + frame. Defaults to (0.0, 0.0, 0.0) + + """ + + def __init__(self, r, phi, z, origin=(0.0, 0.0, 0.0)): + super().__init__() + self.r = r + self.phi = phi + self.z = z + self.origin = origin + + @property + def r(self): + return self._r + + @property + def phi(self): + return self._phi + + @property + def z(self): + return self._z + + @property + def origin(self): + return self._origin + + @r.setter + def r(self, r): + cv.check_type('r coordinate', r, Univariate) + self._r = r + + @phi.setter + def phi(self, phi): + cv.check_type('phi coordinate', phi, Univariate) + self._phi = phi + + @z.setter + def z(self, z): + cv.check_type('z coordinate', z, Univariate) + self._z = z + + @origin.setter + def origin(self, origin): + cv.check_type('origin coordinates', origin, Iterable, Real) + origin = np.asarray(origin) + self._origin = origin + + def to_xml_element(self): + """Return XML representation of the spatial distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing spatial distribution data + + """ + element = ET.Element('space') + element.set('type', 'cylindrical') + element.append(self.r.to_xml_element('r')) + element.append(self.phi.to_xml_element('phi')) + element.append(self.z.to_xml_element('z')) + element.set("origin", ' '.join(map(str, self.origin))) + return element + + @classmethod + def from_xml_element(cls, elem): + """Generate spatial distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.CylindricalIndependent + Spatial distribution generated from XML element + + """ + r = Univariate.from_xml_element(elem.find('r')) + phi = Univariate.from_xml_element(elem.find('phi')) + z = Univariate.from_xml_element(elem.find('z')) + origin = [float(x) for x in elem.get('origin').split()] + return cls(r, phi, z, origin=origin) + class Box(Spatial): """Uniform distribution of coordinates in a rectangular cuboid. diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index a5c45e4e9..12432baf7 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -51,6 +51,72 @@ Position CartesianIndependent::sample(uint64_t* seed) const return {x_->sample(seed), y_->sample(seed), z_->sample(seed)}; } +//============================================================================== +// CylindricalIndependent implementation +//============================================================================== + +CylindricalIndependent::CylindricalIndependent(pugi::xml_node node) +{ + // Read distribution for r-coordinate + if (check_for_node(node, "r")) { + pugi::xml_node node_dist = node.child("r"); + r_ = distribution_from_xml(node_dist); + } else { + // If no distribution was specified, default to a single point at r=0 + double x[] {0.0}; + double p[] {1.0}; + r_ = std::make_unique(x, p, 1); + } + + // Read distribution for phi-coordinate + if (check_for_node(node, "phi")) { + pugi::xml_node node_dist = node.child("phi"); + phi_ = distribution_from_xml(node_dist); + } else { + // If no distribution was specified, default to a single point at phi=0 + double x[] {0.0}; + double p[] {1.0}; + phi_ = std::make_unique(x, p, 1); + } + + // Read distribution for z-coordinate + if (check_for_node(node, "z")) { + pugi::xml_node node_dist = node.child("z"); + z_ = distribution_from_xml(node_dist); + } else { + // If no distribution was specified, default to a single point at z=0 + double x[] {0.0}; + double p[] {1.0}; + z_ = std::make_unique(x, p, 1); + } + + // Read cylinder center coordinates + if (check_for_node(node, "origin")) { + auto origin = get_node_array(node, "origin"); + if (origin.size() == 3) { + origin_ = origin; + } else { + std::stringstream err_msg; + err_msg << "Origin for cylindrical source distribution must be length 3"; + fatal_error(err_msg); + } + } else { + // If no coordinates were specified, default to (0, 0, 0) + origin_ = {0.0, 0.0, 0.0}; + } + +} + +Position CylindricalIndependent::sample(uint64_t* seed) const +{ + double r = r_->sample(seed); + double phi = phi_->sample(seed); + double x = r*cos(phi) + origin_.x; + double y = r*sin(phi) + origin_.y; + double z = z_->sample(seed) + origin_.z; + return {x, y, z}; +} + //============================================================================== // SphericalIndependent implementation //============================================================================== diff --git a/src/source.cpp b/src/source.cpp index 628631fdb..d2e13a173 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -86,6 +86,8 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) type = get_node_value(node_space, "type", true, true); if (type == "cartesian") { space_ = UPtrSpace{new CartesianIndependent(node_space)}; + } else if (type == "cylindrical") { + space_ = UPtrSpace{new CylindricalIndependent(node_space)}; } else if (type == "spherical") { space_ = UPtrSpace{new SphericalIndependent(node_space)}; } else if (type == "box") { diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index 6f05eb3a7..de81a4179 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -63,4 +63,17 @@ 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + + + + + -2.0 0.0 2.0 0.2 0.3 0.2 + + + + + 1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23 + + diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 8f24644f0..eda1a6597 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.004769E-01 3.944044E-03 +2.999424E-01 9.520402E-03 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 5fd06286b..91316acc6 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -38,7 +38,10 @@ class SourceTestHarness(PyAPITestHarness): spatial3 = openmc.stats.Point([1.2, -2.3, 0.781]) spatial4 = openmc.stats.SphericalIndependent(r_dist, theta_dist, phi_dist, - origin=(1.0, 1.0, 0.0)) + origin=(1., 1., 0.)) + spatial5 = openmc.stats.CylindricalIndependent(r_dist, phi_dist, + z_dist, + origin=(1., 1., 0.)) mu_dist = openmc.stats.Discrete([-1., 0., 1.], [0.5, 0.25, 0.25]) phi_dist = openmc.stats.Uniform(0., 6.28318530718) @@ -57,12 +60,13 @@ class SourceTestHarness(PyAPITestHarness): source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3) source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1) source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1) + source5 = openmc.Source(spatial5, angle3, energy3, strength=0.1) settings = openmc.Settings() settings.batches = 10 settings.inactive = 5 settings.particles = 1000 - settings.source = [source1, source2, source3, source4] + settings.source = [source1, source2, source3, source4, source5] settings.export_to_xml()