crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Change logic to run CMFD at cmfd_begin instead of tally_begin

Browse source
This commit is contained in:
Shikhar Kumar 2019-07-22 16:51:29 -04:00
parent bc363cc934
commit 093f90fc64
12 changed files with 120 additions and 202 deletions
Download patch file Download diff file Expand all files Collapse all files

223
openmc/cmfd.py
View file

@ -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

8
tests/regression_tests/cmfd_feed/test.py
View file

@ -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]

2
tests/regression_tests/cmfd_feed_2g/test.py
View file

@ -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

25
tests/regression_tests/cmfd_feed_expanding_window/results_true.dat
View file

@ -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

2
tests/regression_tests/cmfd_feed_expanding_window/test.py
View file

@ -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'

2
tests/regression_tests/cmfd_feed_ng/test.py
View file

@ -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

25
tests/regression_tests/cmfd_feed_ref_d/results_true.dat
View file

@ -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

2
tests/regression_tests/cmfd_feed_ref_d/test.py
View file

@ -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'

25
tests/regression_tests/cmfd_feed_rolling_window/results_true.dat
View file

@ -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

2
tests/regression_tests/cmfd_feed_rolling_window/test.py
View file

@ -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'

2
tests/regression_tests/cmfd_nofeed/test.py
View file

@ -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]

4
tests/regression_tests/cmfd_restart/test.py
View file

@ -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]