Finish implementing compute_xs

This commit is contained in:
Shikhar Kumar 2018-08-20 17:56:22 -04:00
parent 465a9ad838
commit 6798cafde1
16 changed files with 293 additions and 571 deletions

View file

@ -1,508 +0,0 @@
k-combined:
1.169891E+00 6.289489E-03
tally 1:
1.173922E+01
1.385461E+01
2.164076E+01
4.699368E+01
2.906462E+01
8.464937E+01
3.382312E+01
1.147095E+02
3.632006E+01
1.323878E+02
3.655413E+01
1.341064E+02
3.347757E+01
1.124264E+02
2.931336E+01
8.607239E+01
2.182947E+01
4.789565E+01
1.147668E+01
1.325716E+01
tally 2:
2.298190E+01
2.667071E+01
1.600292E+01
1.293670E+01
4.268506E+01
9.161216E+01
3.022909E+01
4.598915E+01
5.680399E+01
1.623879E+02
4.033805E+01
8.196263E+01
6.814742E+01
2.331778E+02
4.851618E+01
1.182330E+02
7.392923E+01
2.740255E+02
5.253586E+01
1.384152E+02
7.332860E+01
2.698608E+02
5.227405E+01
1.371810E+02
6.830172E+01
2.340687E+02
4.867159E+01
1.188724E+02
5.885634E+01
1.736180E+02
4.170434E+01
8.719622E+01
4.371848E+01
9.592893E+01
3.106403E+01
4.844308E+01
2.338413E+01
2.752467E+01
1.636713E+01
1.347770E+01
tally 3:
1.538752E+01
1.196478E+01
1.079685E+00
6.010786E-02
2.911906E+01
4.269070E+01
1.822657E+00
1.671850E-01
3.885421E+01
7.608218E+01
2.541516E+00
3.262451E-01
4.673300E+01
1.097036E+02
2.885307E+00
4.214444E-01
5.059247E+01
1.283984E+02
3.222796E+00
5.237329E-01
5.034856E+01
1.272538E+02
3.230225E+00
5.273424E-01
4.688476E+01
1.103152E+02
2.941287E+00
4.363749E-01
4.013746E+01
8.077506E+01
2.634234E+00
3.520270E-01
2.996887E+01
4.509953E+01
1.946504E+00
1.919104E-01
1.575260E+01
1.248707E+01
1.020705E+00
5.413569E-02
tally 4:
3.049469E+00
4.677325E-01
0.000000E+00
0.000000E+00
2.770358E+00
3.879191E-01
5.514939E+00
1.528899E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.514939E+00
1.528899E+00
2.770358E+00
3.879191E-01
5.032131E+00
1.275040E+00
7.294002E+00
2.675589E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.294002E+00
2.675589E+00
5.032131E+00
1.275040E+00
7.036008E+00
2.490718E+00
8.668860E+00
3.776102E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.668860E+00
3.776102E+00
7.036008E+00
2.490718E+00
8.352414E+00
3.501945E+00
9.345868E+00
4.380719E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.345868E+00
4.380719E+00
8.352414E+00
3.501945E+00
9.093766E+00
4.158282E+00
9.223771E+00
4.270120E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.223771E+00
4.270120E+00
9.093766E+00
4.158282E+00
9.219150E+00
4.264346E+00
8.530966E+00
3.651778E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
8.530966E+00
3.651778E+00
9.219150E+00
4.264346E+00
8.690373E+00
3.785262E+00
7.204424E+00
2.604203E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.204424E+00
2.604203E+00
8.690373E+00
3.785262E+00
7.513640E+00
2.833028E+00
5.326721E+00
1.426975E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.326721E+00
1.426975E+00
7.513640E+00
2.833028E+00
5.662215E+00
1.607757E+00
2.848381E+00
4.093396E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.848381E+00
4.093396E-01
5.662215E+00
1.607757E+00
3.025812E+00
4.597241E-01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
tally 5:
1.538652E+01
1.196332E+01
2.252427E+00
2.605738E-01
2.911344E+01
4.267319E+01
3.873926E+00
7.615035E-01
3.884516E+01
7.604619E+01
5.280610E+00
1.414008E+00
4.672391E+01
1.096625E+02
6.261805E+00
1.983205E+00
5.058447E+01
1.283588E+02
6.733810E+00
2.278242E+00
5.033589E+01
1.271898E+02
6.714658E+00
2.273652E+00
4.687563E+01
1.102719E+02
6.215002E+00
1.956978E+00
4.013134E+01
8.075062E+01
5.253064E+00
1.396224E+00
2.996497E+01
4.508840E+01
3.818076E+00
7.509442E-01
1.574994E+01
1.248291E+01
2.219928E+00
2.515492E-01
cmfd indices
1.000000E+01
1.000000E+00
1.000000E+00
1.000000E+00
k cmfd
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.170416E+00
1.172966E+00
1.165537E+00
1.170979E+00
1.161922E+00
1.157523E+00
1.158873E+00
1.162877E+00
1.167101E+00
1.168130E+00
1.170570E+00
1.168115E+00
1.174081E+00
1.169458E+00
1.167848E+00
1.165116E+00
cmfd entropy
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.203643E+00
3.207943E+00
3.213367E+00
3.214360E+00
3.219634E+00
3.222232E+00
3.221744E+00
3.224544E+00
3.225990E+00
3.227769E+00
3.227417E+00
3.230728E+00
3.231662E+00
3.233316E+00
3.233193E+00
3.232564E+00
cmfd balance
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
4.009062E-03
4.431773E-03
3.152666E-03
3.510383E-03
2.052089E-03
2.068651E-03
1.502427E-03
1.589825E-03
1.566020E-03
1.219160E-03
1.017888E-03
9.771622E-04
1.010120E-03
1.073382E-03
1.172758E-03
9.827332E-04
cmfd dominance ratio
0.000E+00
0.000E+00
0.000E+00
0.000E+00
5.397E-01
5.425E-01
5.481E-01
5.473E-01
5.503E-01
5.502E-01
5.483E-01
5.520E-01
5.505E-01
3.216E-01
5.373E-01
5.517E-01
5.508E-01
5.524E-01
5.524E-01
5.523E-01
cmfd openmc source comparison
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
6.959834E-03
5.655657E-03
3.886185E-03
4.035116E-03
3.043277E-03
5.455475E-03
4.515311E-03
2.439840E-03
2.114032E-03
2.673132E-03
2.431749E-03
4.330928E-03
3.404647E-03
3.680298E-03
3.309620E-03
3.705541E-03
cmfd source
4.697085E-02
7.920706E-02
1.107968E-01
1.250932E-01
1.383930E-01
1.380648E-01
1.246874E-01
1.113705E-01
8.203754E-02
4.337882E-02

View file

@ -7,6 +7,8 @@ cmfd_mesh.lower_left = [-10, -1, -1]
cmfd_mesh.upper_right = [10, 1, 1]
cmfd_mesh.dimension = [10, 1, 1]
cmfd_mesh.albedo = [0., 0., 1., 1., 1., 1.]
cmfd_mesh.energy = [0., 10000., 10000000.]
cmfd_mesh.map = [0,1,1,1,1,1,1,1,1,0]
# Initialize CMFDRun object
cmfd_run = openmc.CMFDRun()
@ -14,7 +16,10 @@ cmfd_run = openmc.CMFDRun()
# Set all runtime parameters (cmfd_mesh, tolerances, tally_resets, etc)
# All error checking done under the hood when setter function called
cmfd_run.cmfd_mesh = cmfd_mesh
cmfd_run.cmfd_reset = [5,10]
cmfd_run.cmfd_begin = 5
cmfd_run.cmfd_display = 'dominance'
cmfd_run.cmfd_feedback = True
cmfd_run.gauss_seidel_tolerance = [1.e-15, 1.e-20]
# Run CMFD
cmfd_run.run()

View file

@ -6,10 +6,11 @@
<upper_right> 10 1 1 </upper_right>
<dimension> 10 1 1 </dimension>
<albedo> 0.0 0.0 1.0 1.0 1.0 1.0 </albedo>
<energy> 0.0 10000.0 10000000.0</energy>
<map> 1 2 2 2 2 2 2 2 2 1</map>
</mesh>
<tally_reset> 5 10 </tally_reset>
<begin>5</begin>
<display>dominance</display>
<feedback>true</feedback>
<feedback>false</feedback>
<gauss_seidel_tolerance>1.e-15 1.e-20</gauss_seidel_tolerance>
</cmfd>

View file

@ -1,4 +1,5 @@
from ctypes import c_int, c_int32, c_int64, c_double, c_char_p, POINTER
from ctypes import (c_int, c_int32, c_int64, c_double, c_char_p, c_bool,
POINTER)
from . import _dll
from .core import _DLLGlobal
@ -17,9 +18,13 @@ _dll.openmc_get_seed.restype = c_int64
class _Settings(object):
# Attributes that are accessed through a descriptor
batches = _DLLGlobal(c_int32, 'n_batches')
current_batch = _DLLGlobal(c_int, 'openmc_current_batch')
generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch')
inactive = _DLLGlobal(c_int32, 'n_inactive')
master = _DLLGlobal(c_bool, 'openmc_master')
particles = _DLLGlobal(c_int64, 'n_particles')
restart_batch = _DLLGlobal(c_int, 'openmc_restart_batch')
restart_run = _DLLGlobal(c_bool, 'openmc_restart_run')
verbosity = _DLLGlobal(c_int, 'openmc_verbosity')
@property

View file

@ -10,6 +10,7 @@ References
"""
# TODO: Check to make sure no redundant import statements
from collections.abc import Iterable
from numbers import Real, Integral
from xml.etree import ElementTree as ET
@ -21,9 +22,28 @@ from openmc.clean_xml import clean_xml_indentation
from openmc.checkvalue import (check_type, check_length, check_value,
check_greater_than, check_less_than)
# Maximum/minimum neutron energies, from src/api.F90
ENERGY_MAX_NEUTRON = np.inf
ENERGY_MIN_NEUTRON = 0.
"""
--------------
CMFD CONSTANTS
--------------
"""
# Maximum/minimum neutron energies
_ENERGY_MAX_NEUTRON = np.inf
_ENERGY_MIN_NEUTRON = 0.
# Tolerance for detecting zero flux values
_TINY_BIT = 1.e-8
# For non-accelerated regions on coarse mesh overlay
_CMFD_NOACCEL = 99999
# Constant to represent a zero flux "albedo"
_ZERO_FLUX = 999.0
# Constant for writing out no residual
_CMFD_NORES = 99999.0
class CMFDMesh(object):
"""A structured Cartesian mesh used for Coarse Mesh Finite Difference (CMFD)
@ -156,7 +176,7 @@ class CMFDMesh(object):
def map(self, meshmap):
check_type('CMFD mesh map', meshmap, Iterable, Integral)
for m in meshmap:
check_value('CMFD mesh map', m, [1, 2])
check_value('CMFD mesh map', m, [0, 1])
self._map = meshmap
def _get_xml_element(self):
@ -186,7 +206,7 @@ class CMFDMesh(object):
if self.map is not None:
subelement = ET.SubElement(element, "map")
subelement.text = ' '.join(map(str, self.map))
subelement.text = ' '.join(map(str, [self.map[i]+1 for i in range(len(self.map))]))
return element
@ -409,8 +429,8 @@ class CMFD(object):
# Check upper right coordinates are greater than lower left
if np.any(np.array(mesh.upper_right) <= np.array(mesh.lower_left)):
raise ValueError('CMFD mesh requires upper right '
'coordinates to be greater than lower '
'left coordinates')
'coordinates to be greater than lower '
'left coordinates')
mesh.width = np.true_divide(
(np.array(mesh.upper_right) - np.array(mesh.lower_left)), \
np.array(mesh.dimension))
@ -577,11 +597,9 @@ class CMFDRun(object):
Attributes
----------
To add: cmfd_coremap: Flag for active core map
indices: Stores spatial and group dimensions as [nx, ny, nz, ng]
To add: indices: Stores spatial and group dimensions as [nx, ny, nz, ng]
egrid: energy grid used for CMFD acceleration
albedo: Albedo for global boundary conditions, taken from CMFD mesh. Set to [1,1,1,1,1,1] if not specified by user
cmfd_coremap: Optional acceleration map to overlay on coarse mesh spatial grid, taken from CMFD mesh
n_cmfd_resets: Number of elements in tally_reset, list that stores batches where CMFD tallies should be reset
cmfd_atoli: Absolute GS tolerance, set by gauss_seidel_tolerance
cmfd_rtoli: Relative GS tolerance, set by gauss_seidel_tolerance
@ -592,7 +610,6 @@ class CMFDRun(object):
Set to true if user specifies energy grid in CMFDMesh, false otherwise
cmfd_on: Boolean to tell if cmfd solver should be initiated, based on whether current batch has reached variable
cmfd_begin
batch_num: Current batch
# Look at cmfd_header.F90 for description
flux
totalxs
@ -613,9 +630,12 @@ class CMFDRun(object):
src_cmp
dom
k_cmfd
keff_bal
mat_dim
TODO: Put descriptions for all methods in CMFDRun
TODO All timing variables
TODO Get rid of CMFD constants
"""
@ -641,7 +661,6 @@ class CMFDRun(object):
self._cmfd_write_matrices = False
# External variables used during runtime but users don't have control over
self._cmfd_coremap = False
self._indices = np.zeros(4, dtype=int)
self._egrid = None
self._albedo = None
@ -654,7 +673,8 @@ class CMFDRun(object):
self._cmfd_tally_ids = None
self._energy_filters = None
self._cmfd_on = False
self._batch_num = 1
self._mat_dim = _CMFD_NOACCEL
self._keff_bal = None
# Numpy arrays used to build CMFD matrices
self._flux = None
@ -745,7 +765,7 @@ class CMFDRun(object):
def cmfd_begin(self, cmfd_begin):
check_type('CMFD begin batch', cmfd_begin, Integral)
check_greater_than('CMFD begin batch', cmfd_begin, 0)
self._cmfd_begin = begin
self._cmfd_begin = cmfd_begin
@dhat_reset.setter
def dhat_reset(self, dhat_reset):
@ -793,12 +813,12 @@ class CMFDRun(object):
# Check lower left defined
if mesh.lower_left is None:
raise ValueError('CMFD mesh requires lower left coordinates '
'to be specified')
'to be specified')
# Check that both upper right and width both not defined
if mesh.upper_right is not None and mesh.width is not None:
raise ValueError('Both upper right coordinates and width '
'cannot be specified for CMFD mesh')
'cannot be specified for CMFD mesh')
# Check that at least one of width or upper right is defined
if mesh.upper_right is None and mesh.width is None:
@ -865,7 +885,8 @@ class CMFDRun(object):
check_type('CMFD write matrices', cmfd_write_matrices, bool)
self._cmfd_write_matrices = cmfd_write_matrices
def run(self):
def run(self, mpi_procs=None, omp_threads=None):
# TODO: Add logic for mpi_procs, omp_threads as args
openmc.capi.init()
self._configure_cmfd()
openmc.capi.simulation_init()
@ -877,13 +898,13 @@ class CMFDRun(object):
# CMFD update or skip it entirely if it is a restart run
status = openmc.capi.next_batch_between_cmfd_init_execute()
if status != 0:
self._cmfd_execute()
if self._cmfd_on:
self._execute_cmfd()
# Status now determines whether another batch should be run
# or simulation should be terminated.
status = openmc.capi.next_batch_after_cmfd_execute()
if status != 0:
break
self._batch_num += 1
openmc.capi.simulation_finalize()
openmc.capi.finalize()
@ -902,12 +923,14 @@ class CMFDRun(object):
self._allocate_cmfd()
def _read_cmfd_input(self):
# TODO: Print message with verbosity
# Print message
print(' Configuring CMFD parameters for simulation')
# Check if CMFD mesh is defined
if self._cmfd_mesh is None:
raise ValueError('No CMFD mesh has been specified for '
'simulation')
'simulation')
# Set spatial dimensions of CMFD object
# Iterate through each element of self._cmfd_mesh.dimension
@ -923,7 +946,7 @@ class CMFDRun(object):
self._energy_filters = True
# TODO: MG mode check
else:
self._egrid = np.array([ENERGY_MIN_NEUTRON, ENERGY_MAX_NEUTRON])
self._egrid = np.array([_ENERGY_MIN_NEUTRON, _ENERGY_MAX_NEUTRON])
self._indices[3] = 1
self._energy_filters = False
@ -933,14 +956,16 @@ class CMFDRun(object):
else:
self._albedo = np.array([1.,1.,1.,1.,1.,1.])
# Get acceleration map
# Get acceleration map, otherwise set all regions to be accelerated
if self._cmfd_mesh.map is not None:
check_length('CMFD coremap', self._cmfd_mesh.map,
np.product(self._indices[0:3]))
self._coremap = np.array(self._cmfd_mesh.map).reshape(( \
self._indices[0], self._indices[1], \
self._indices[2]))
self._cmfd_coremap = True
else:
self._coremap = np.ones((self._indices[0], self._indices[1], \
self._indices[2]))
# Set number of batches where cmfd tallies should be reset
if self._cmfd_reset is not None:
@ -954,9 +979,6 @@ class CMFDRun(object):
self._create_cmfd_tally()
def _allocate_cmfd(self):
# Determine number of batches through C API
n_batches = openmc.capi.settings.batches
# Extract spatial and energy indices
nx = self._indices[0]
ny = self._indices[1]
@ -964,57 +986,92 @@ class CMFDRun(object):
ng = self._indices[3]
# Allocate flux, cross sections and diffusion coefficient
self._flux = np.zeros((ng, nx, ny, nz))
self._totalxs = np.zeros((ng, nx, ny, nz))
self._p1scattxs = np.zeros((ng, nx, ny, nz))
self._scattxs = np.zeros((ng, ng, nx, ny, nz))
self._nfissxs = np.zeros((ng, ng, nx, ny, nz))
self._diffcof = np.zeros((ng, nx, ny, nz))
self._flux = np.zeros((nx, ny, nz, ng))
self._totalxs = np.zeros((nx, ny, nz, ng))
self._p1scattxs = np.zeros((nx, ny, nz, ng))
self._scattxs = np.zeros((nx, ny, nz, ng, ng)) # Outgoing, incoming
self._nfissxs = np.zeros((nx, ny, nz, ng, ng)) # Outgoing, incoming
self._diffcof = np.zeros((nx, ny, nz, ng))
# Allocate dtilde and dhat
self._dtilde = np.zeros((6, ng, nx, ny, nz))
self._dhat = np.zeros((6, ng, nx, ny, nz))
self._dtilde = np.zeros((6, nx, ny, nz, ng))
self._dhat = np.zeros((6, nx, ny, nz, ng))
# Allocate dimensions for each box (here for general case)
self._hxyz = np.zeros((3, nx, ny, nz))
# Allocate surface currents
self._current = np.zeros((12, ng, nx, ny, nz))
self._current = np.zeros((nx, ny, nz, ng, 12))
# Allocate source distributions
self._cmfd_src = np.zeros((ng, nx, ny, nz))
self._openmc_src = np.zeros((ng, nx, ny, nz))
self._cmfd_src = np.zeros((nx, ny, nz, ng))
self._openmc_src = np.zeros((nx, ny, nz, ng))
# Allocate source weight modification variables
self._sourcecounts = np.zeros((ng, nx*ny*nz))
self._weightfactors = np.ones((ng, nx, ny, nz))
self._sourcecounts = np.zeros((nx*ny*nz, ng))
self._weightfactors = np.ones((nx, ny, nz, ng))
# Allocate batchwise parameters
self._entropy = np.zeros((n_batches, ))
self._balance = np.zeros((n_batches, ))
self._src_cmp = np.zeros((n_batches, ))
self._dom = np.zeros((n_batches, ))
self._k_cmfd = np.zeros((n_batches, ))
self._entropy = []
self._balance = []
self._src_cmp = []
self._dom = []
self._k_cmfd = []
def _cmfd_init_batch(self):
if self._cmfd_begin == self._batch_num:
current_batch = openmc.capi.settings.current_batch
restart_run = openmc.capi.settings.restart_run
restart_batch = openmc.capi.settings.restart_batch
if self._cmfd_begin == current_batch:
self._cmfd_on = True
# TODO: Figure out how to tell if part of restart run
# if restart_run:
# return
# If this is a restart run and we are just replaying batches leave
# if (restart_run .and. current_batch <= restart_batch) return
# TODO: Test restart_batch
if restart_run and current_batch <= restart_batch:
return
# Check to reset tallies
if self._n_cmfd_resets > 0 and self._batch_num in self._cmfd_reset:
if (self._n_cmfd_resets > 0
and current_batch in self._cmfd_reset):
self._cmfd_tally_reset()
def _cmfd_execute(self):
pass
def _execute_cmfd(self):
# CMFD single processor on master
if openmc.capi.settings.master:
# TODO
#! Start cmfd timer
#call time_cmfd % start()
# Create cmfd data from OpenMC tallies
self._set_up_cmfd()
'''
! Call solver
call cmfd_solver_execute()
! Save k-effective
cmfd % k_cmfd(current_batch) = cmfd % keff
! check to perform adjoint on last batch
if (current_batch == n_batches .and. cmfd_run_adjoint) then
call cmfd_solver_execute(adjoint=.true.)
end if
end if
! calculate fission source
call calc_fission_source()
! calculate weight factors
call cmfd_reweight(.true.)
! stop cmfd timer
if (master) call time_cmfd % stop()
'''
def _cmfd_tally_reset(self):
# TODO: How does write_message in error.F90 work?
# TODO: Print message with verbosity
# Print message
print(' CMFD tallies reset')
@ -1023,6 +1080,168 @@ class CMFDRun(object):
for tally_id in self._cmfd_tally_ids:
tallies[tally_id].reset()
def _set_up_cmfd(self):
# Check for core map and set it up
if (self._mat_dim == _CMFD_NOACCEL):
self._set_coremap()
# Calculate all cross sections based on reaction rates from last batch
self._compute_xs()
'''
! Compute effective downscatter cross section
if (cmfd_downscatter) call compute_effective_downscatter()
! Check neutron balance
call neutron_balance()
! Calculate dtilde
call compute_dtilde()
! Calculate dhat
call compute_dhat()
'''
def _set_coremap(self):
self._mat_dim = np.sum(self._coremap)
self._coremap = np.where(self._coremap == 0,
_CMFD_NOACCEL, self._coremap)
def _compute_xs(self):
# Extract energy indices
ng = self._indices[3]
# Set flux object and source distribution all to zeros
self._flux.fill(0.)
self._openmc_src.fill(0.)
# Reset keff_bal to zero
self._keff_bal = 0.
# Get tallies in-memory
tallies = openmc.capi.tallies
# Set conditional numpy array as boolean vector based on coremap
# Repeat each value for number of groups in problem
is_cmfd_accel = np.repeat(self._coremap.ravel() != _CMFD_NOACCEL, ng)
# Get flux from CMFD tally 0
tally_id = self._cmfd_tally_ids[0]
tally_results = tallies[tally_id].results[:,0,1]
flux = np.where(is_cmfd_accel, tally_results, 0.)
# Detect zero flux, abort if located
if np.any(flux[is_cmfd_accel] < _TINY_BIT):
# Get index of zero flux in flux array
idx = np.argmax(np.where(is_cmfd_accel, flux, 1) < _TINY_BIT)
# Convert scalar idx to index in flux matrix
mat_idx = np.unravel_index(idx, self._flux.shape)
# Throw error message (one-based indexing)
# Index of group is flipped
err_message = 'Detected zero flux without coremap overlay' + \
' at mesh: (' + \
', '.join(str(i+1) for i in mat_idx[:-1]) + \
') in group ' + str(ng-mat_idx[-1])
raise ValueError(err_message)
# Store flux and reshape
# Flux is flipped in energy axis as tally results are given in reverse
# order of energy group
self._flux = np.flip(flux.reshape(self._flux.shape), axis=3)
# Get total rr and convert to total xs from CMFD tally 0
tally_results = tallies[tally_id].results[:,1,1]
totalxs = np.divide(tally_results, flux, \
where=flux>0, out=np.zeros_like(tally_results))
# Store total xs and reshape
# Total xs is flipped in energy axis as tally results are given in
# reverse order of energy group
self._totalxs = np.flip(totalxs.reshape(self._totalxs.shape), axis=3)
# Get scattering xs from CMFD tally 1
# flux is repeated to account for extra dimensionality of scattering xs
tally_id = self._cmfd_tally_ids[1]
tally_results = tallies[tally_id].results[:,0,1]
scattxs = np.divide(tally_results, \
np.repeat(flux, ng), \
where=np.repeat(flux>0, ng), \
out=np.zeros_like(tally_results))
# Store scattxs and reshape
# Scattering xs is flipped in both incoming and outgoing energy axes
# as tally results are given in reverse order of energy group
self._scattxs = np.flip(scattxs.reshape(self._scattxs.shape), axis=3)
self._scattxs = np.flip(self._scattxs.reshape(self._scattxs.shape), \
axis=4)
# Get nu-fission xs from CMFD tally 1
# flux is repeated to account for extra dimensionality of nu-fission xs
tally_results = tallies[tally_id].results[:,1,1]
num_realizations = tallies[tally_id].num_realizations
nfissxs = np.divide(tally_results, \
np.repeat(flux, ng), \
where=np.repeat(flux>0, ng), \
out=np.zeros_like(tally_results))
# Store nfissxs and reshape
# Nu-fission xs is flipped in both incoming and outgoing energy axes
# as tally results are given in reverse order of energy group
self._nfissxs = np.flip(nfissxs.reshape(self._nfissxs.shape), axis=3)
self._nfissxs = np.flip(self._nfissxs.reshape(self._nfissxs.shape), \
axis=4)
# Filter nu-fission tally results to compute openmc source distribution
tally_results = np.where(np.repeat(flux>0, ng), tally_results, \
0.)
# Openmc source distribution is sum of nu-fission rr in outgoing energies
openmc_src = np.sum(tally_results.reshape(self._nfissxs.shape),
axis=3)
# Store openmc_src
# Openmc source is flipped in energy axis as tally results are given
# in reverse order of energy group
self._openmc_src = np.flip(openmc_src, axis=3)
# Compute k_eff from source distribution
self._keff_bal = np.sum(self._openmc_src) / num_realizations
# Normalize openmc source distribution
self._openmc_src /= np.sum(self._openmc_src) * self._norm
# Get surface currents from CMFD tally 2
tally_id = self._cmfd_tally_ids[2]
tally_results = tallies[tally_id].results[:,0,1]
# Filter tally results to include only accelerated regions
tally_results = np.where(np.repeat(flux>0, 12), tally_results, 0.)
# Reshape and store current
# Current is flipped in energy axis as tally results are given in
# reverse order of energy group
self._current = np.flip(tally_results.reshape(self._current.shape), \
axis=3)
# Get p1 scatter xs from CMFD tally 3
tally_id = self._cmfd_tally_ids[3]
tally_results = tallies[tally_id].results[:,0,1]
# Reshape and extract only p1 data from tally results (no need for p0 data)
p1scattrr = tally_results.reshape(self._p1scattxs.shape+(2,))[:,:,:,1]
# Store p1 scatter xs
# p1 scatter xs is flipped in energy axis as tally results are given in
# reverse order of energy group
self._p1scattxs = np.divide(np.flip(p1scattrr, axis=3), self._flux, \
where=self._flux>0, \
out=np.zeros_like(p1scattrr))
# Calculate and store diffusion coefficient
self._diffcof = np.where(self._flux > 0, 1.0 / (3.0 * \
(self._totalxs - self._p1scattxs)), 0.)
def _create_cmfd_tally(self):
# Create Mesh object based on CMFDMesh, stored internally
cmfd_mesh = openmc.capi.Mesh()
@ -1083,7 +1302,7 @@ class CMFDRun(object):
tally.filters = [mesh_filter, energy_filter]
else:
tally.filters = [mesh_filter]
# Set scores for tally
# Set scores, type, and estimator for tally
tally.scores = ['flux', 'total']
tally.type = 'volume'
tally.estimator = 'analog'
@ -1095,7 +1314,7 @@ class CMFDRun(object):
tally.filters = [mesh_filter, energy_filter, energyout_filter]
else:
tally.filters = [mesh_filter]
# Set scores for tally
# Set scores, type, and estimator for tally
tally.scores = ['nu-scatter', 'nu-fission']
tally.type = 'volume'
tally.estimator = 'analog'
@ -1107,7 +1326,7 @@ class CMFDRun(object):
tally.filters = [meshsurface_filter, energy_filter]
else:
tally.filters = [meshsurface_filter]
# Set scores for tally
# Set scores, type, and estimator for tally
tally.scores = ['current']
tally.type = 'mesh-surface'
tally.estimator = 'analog'

View file

@ -254,7 +254,7 @@ contains
! Set the energy bin if needed
if (energy_filters) then
filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
filter_matches(i_filter_ein) % bins % data(1) = 12*(ng - h) + 1
end if
score_index = 0

View file

@ -65,7 +65,7 @@ module simulation_header
! ============================================================================
! MISCELLANEOUS VARIABLES
integer :: restart_batch
integer(C_INT), bind(C, name='openmc_restart_batch') :: restart_batch
! Flag for enabling cell overlap checking during transport
integer(8), allocatable :: overlap_check_cnt(:)