Clean up comments, add vector_write function in Fortran, create vectorized calc_fission_source

This commit is contained in:
Shikhar Kumar 2018-08-26 22:45:58 -04:00
parent 961178ed80
commit c2858bc941
3 changed files with 194 additions and 84 deletions

View file

@ -16,6 +16,7 @@ from xml.etree import ElementTree as ET # TODO Remove
import sys # TODO Remove
import numpy as np
from scipy import sparse
import time
import openmc.capi
from openmc.clean_xml import clean_xml_indentation # TODO Remove
@ -637,11 +638,9 @@ class CMFDRun(object):
mat_dim
TODO: Put descriptions for all methods in CMFDRun
TODO All timing variables
TODO Get rid of CMFD constants
TODO Get rid of unused variables defined in init
TODO Make sure all self variables defined in init
TODO Clean up logic for adjoint, understand what different adjoint types are doing
TODO Check to make sure no compatibility issues with numpy arrays for input variables
TODO Create write_vector function in cmfd_solver.F90
@ -679,7 +678,7 @@ class CMFDRun(object):
self._cmfd_on = False
self._mat_dim = _CMFD_NOACCEL
self._keff_bal = None
self._cmfd_adjoint_type = None
self._cmfd_adjoint_type = "physical"
self._keff = None
self._adj_keff = None
self._phi = None
@ -708,6 +707,10 @@ class CMFDRun(object):
self._k_cmfd = None
self._resnb = None
self._time_cmfd = None
self._time_cmfdbuild = None
self._time_cmfdsolve = None
@property
def cmfd_begin(self):
return self._cmfd_begin
@ -740,6 +743,10 @@ class CMFDRun(object):
def norm(self):
return self._norm
@property
def cmfd_adjoint_type(self):
return self._cmfd_adjoint_type
@property
def cmfd_power_monitor(self):
return self._cmfd_power_monitor
@ -846,6 +853,13 @@ class CMFDRun(object):
check_type('CMFD norm', norm, Real)
self._norm = norm
@cmfd_adjoint_type.setter
def cmfd_adjoint_type(self, adjoint_type):
check_type('CMFD adjoint type', adjoint_type, str)
check_value('CMFD adjoint type', adjoint_type,
['math', 'phyical'])
self._cmfd_adjoint_type = adjoint_type
@cmfd_power_monitor.setter
def cmfd_power_monitor(self, cmfd_power_monitor):
check_type('CMFD power monitor', cmfd_power_monitor, bool)
@ -877,47 +891,65 @@ class CMFDRun(object):
self._cmfd_write_matrices = cmfd_write_matrices
def run(self, mpi_procs=None, omp_num_threads=None):
# Check number of OpenMP threads is valid input and initialize C API
if omp_num_threads is not None:
check_type('OpenMP num threads', omp_num_threads, Integral)
openmc.capi.init(args=['-s',str(omp_num_threads)])
else:
openmc.capi.init()
# Configure cmfd parameters and tallies
self._configure_cmfd()
# Initialize simulation
openmc.capi.simulation_init()
while(True):
# Run everything in next batch before initializing cmfd
openmc.capi.next_batch_before_cmfd_init()
# Initialize CMFD batch
self._cmfd_init_batch()
# Run everything in next batch in between initializing and
# executing CMFD
status = openmc.capi.next_batch_between_cmfd_init_execute()
# Status determines whether batch should continue with a
# CMFD update or skip it entirely if it is a restart run
status = openmc.capi.next_batch_between_cmfd_init_execute()
if status != 0:
# Perform CMFD calculation if on
if self._cmfd_on:
self._execute_cmfd()
# Status now determines whether another batch should be run
# or simulation should be terminated.
# Run everything in next batch after executing 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
# Finalize simuation
openmc.capi.simulation_finalize()
# Finalize and free memory
openmc.capi.finalize()
def _configure_cmfd(self):
# Read in cmfd input from python
# Read in cmfd input defined in Python
self._read_cmfd_input()
# TODO
# Initialize timers
#call time_cmfd % reset()
#call time_cmfdbuild % reset()
#call time_cmfdsolve % reset()
self._time_cmfd = 0.0
self._time_cmfdbuild = 0.0
self._time_cmfdsolve = 0.0
# Initialize all numpy arrays used for cmfd solver
self._allocate_cmfd()
def _read_cmfd_input(self):
# Print message
# Print message to user
if openmc.capi.settings.verbosity >= 7 and openmc.capi.settings.master:
print(' Configuring CMFD parameters for simulation')
@ -985,6 +1017,7 @@ class CMFDRun(object):
self._dhat = np.zeros((nx, ny, nz, ng, 6))
# Allocate dimensions for each box (assume fixed mesh dimensions)
# TODO Update this
self._hxyz = np.zeros((3))
# Allocate surface currents
@ -1006,14 +1039,18 @@ class CMFDRun(object):
self._k_cmfd = []
def _cmfd_init_batch(self):
# Get simulation parameters through C API
current_batch = openmc.capi.settings.current_batch
restart_run = openmc.capi.settings.restart_run
restart_batch = openmc.capi.settings.restart_batch
# Check to activate CMFD diffusion and possible feedback
if self._cmfd_begin == current_batch:
self._cmfd_on = True
# TODO: Test restart_batch
# If this is a restart run we are just replaying batches so don't
# execute anything
if restart_run and current_batch <= restart_batch:
return
@ -1023,11 +1060,10 @@ class CMFDRun(object):
self._cmfd_tally_reset()
def _execute_cmfd(self):
# CMFD single processor on master
# Run CMFD on single processor on master
if openmc.capi.settings.master:
# TODO
#! Start cmfd timer
#call time_cmfd % start()
time_start_cmfd = time.time()
# Create cmfd data from OpenMC tallies
self._set_up_cmfd()
@ -1037,22 +1073,22 @@ class CMFDRun(object):
# Save k-effective
self._k_cmfd.append(self._keff)
'''
! TODO 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
'''
# calculate fission source
# Check to perform adjoint on last batch
if (openmc.capi.settings.current_batch == \
openmc.capi.settings.batches and self._cmfd_run_adjoint):
self._cmfd_solver_execute(adjoint=True)
# Calculate fission source
self._calc_fission_source()
'''
! TODO calculate weight factors
! TODO Calculate weight factors
call cmfd_reweight(.true.)
! TODO stop cmfd timer
if (master) call time_cmfd % stop()
'''
# Stop cmfd timer
if openmc.capi.settings.master:
time_stop_cmfd = time.time()
self._time_cmfd += time_stop_cmfd - time_start_cmfd
def _cmfd_tally_reset(self):
# Print message
@ -1065,9 +1101,8 @@ class CMFDRun(object):
tallies[tally_id].reset()
def _set_up_cmfd(self):
# Check for core map and set it up
# Set up CMFD coremap
if (self._mat_dim == _CMFD_NOACCEL):
# TODO: Don't reshape coremap before this point, reshape in set_coremap
self._set_coremap()
# Calculate all cross sections based on reaction rates from last batch
@ -1087,30 +1122,31 @@ class CMFDRun(object):
# Calculate dhat
self._compute_dhat()
# Calculate dhat
self._compute_dhat2()
def _cmfd_solver_execute(self, adjoint=False):
# TODO Check for physical adjoint
# Check for physical adjoint
physical_adjoint = adjoint and self._cmfd_adjoint_type == 'physical'
# TODO Start timer for build
# call time_cmfdbuild % start()
# Start timer for build
time_start_buildcmfd = time.time()
# Initialize matrices and vectors
# Build loss and production matrices
loss, prod = self._build_matrices(physical_adjoint)
# TODO Check for mathematical adjoint calculation
#if adjoint_calc and self._cmfd_adjoint_type == 'math':
# self._compute_adjoint()
# Check for mathematical adjoint calculation
if adjoint and self._cmfd_adjoint_type == 'math':
loss, prod = self._compute_adjoint(loss, prod)
# TODO Stop timer for build
# call time_cmfdbuild % stop()
# Stop timer for build
time_stop_buildcmfd = time.time()
self._time_cmfdbuild += time_stop_buildcmfd - time_start_buildcmfd
# Begin power iteration
# TODO call time_cmfdsolve % start()
time_start_solvecmfd = time.time()
phi, keff, dom = self._execute_power_iter(loss, prod)
# TODO call time_cmfdsolve % stop()
time_stop_solvecmfd = time.time()
self._time_cmfdsolve += time_stop_solvecmfd - time_start_solvecmfd
# Save results, normalizing phi to sum to 1
if adjoint:
@ -1135,20 +1171,61 @@ class CMFDRun(object):
'''
def _calc_fission_source(self):
pass
# Extract number of groups and number of accelerated regions
nx = self._indices[0]
ny = self._indices[1]
nz = self._indices[2]
ng = self._indices[3]
n = self._mat_dim
# Reset cmfd source to 0
self._cmfd_src.fill(0.)
# Calculate volume
vol = np.product(self._hxyz)
# Reshape phi by number of groups
phi = self._phi.reshape((n, ng))
# Extract indices of coremap that are accelerated
idx = np.where(self._coremap != _CMFD_NOACCEL)
# Initialize CMFD flux map that maps phi to actualy spatial and group
# indices of problem
cmfd_flux = np.zeros((nx, ny, nz, ng))
# Loop over all groups and set CMFD flux based on indices of coremap
# and values of phi
for g in range(ng):
flux[idx + (np.full((n,),g),)] = phi[:,g]
# Compute fission source
self._cmfd_src = np.sum(self._nfissxs[:,:,:,:,:] * \
cmfd_flux[:,:,:,:,np.newaxis], axis=3) * vol
def _build_matrices(self, adjoint):
# Set up matrices
# Build loss matrix
loss = self._build_loss_matrix(adjoint)
# Build production matrix
prod = self._build_prod_matrix(adjoint)
'''
# TODO Write matrices
if (cmfd_write_matrices) then
call loss % write('loss.dat')
call prod % write('prod.dat')
end if
'''
# TODO Write out matrices
#if self._cmfd_write_matrices:
# self._write_matrix(loss, 'loss.dat')
# self._write_matrix(prod, 'prod.dat')
return loss, prod
def _compute_adjoint(self, loss, prod):
# Transpose matrices
loss = np.transpose(loss)
prod = np.transpose(prod)
# TODO Write out matrices
#if self._cmfd_write_matrices:
# self._write_matrix(loss, 'adj_loss.dat')
# self._write_matrix(prod, 'adj_prod.dat')
return loss, prod
@ -1170,6 +1247,7 @@ class CMFDRun(object):
jo = np.zeros((6,))
for irow in range(n):
# Get indices for row in matrix
i,j,k,g = self._matrix_to_indices(irow, nx, ny, nz, ng)
# Retrieve cell data
@ -1219,6 +1297,7 @@ class CMFDRun(object):
val = jnet + totxs - scattxsgg
loss[irow, irow] = val
# Begin loop over off diagonal in-scattering
for h in range(ng):
# Cycle though if h=g, value already banked in removal xs
if h == g:
@ -1227,14 +1306,13 @@ class CMFDRun(object):
# Get neighbor matrix index
scatt_mat_idx = self._indices_to_matrix(i,j,k, h, ng)
# TODO Check for adjoint
#if (adjoint_calc) then
#! Get scattering macro xs, transposed!
#scattxshg = cmfd%scattxs(g, h, i, j, k)
#else
# Get scattering macro xs
scattxshg = self._scattxs[i, j, k, h, g]
#end if
# Check for adjoint
if adjoint:
# Get scattering macro xs, transposed!
scattxshg = self._scattxs[i, j, k, g, h]
else:
# Get scattering macro xs
scattxshg = self._scattxs[i, j, k, h, g]
# Negate the scattering xs
val = -1.0*scattxshg
@ -1268,24 +1346,26 @@ class CMFDRun(object):
prod = np.zeros((n, n))
for irow in range(n):
# Get indices for row in matrix
i,j,k,g = self._matrix_to_indices(irow, nx, ny, nz, ng)
# Check if at a reflector
if self._coremap[i,j,k] == _CMFD_NOACCEL:
continue
# loop around all other groups
# Loop around all other groups
for h in range(ng):
# Get matrix column location
hmat_idx = self._indices_to_matrix(i,j,k, h, ng)
# TODO check for adjoint and bank val
#if (adjoint_calc) then
# ! get nu-fission cross section from cell
# nfissxs = cmfd%nfissxs(g,h,i,j,k)
#else
# get nu-fission cross section from cell
nfissxs = self._nfissxs[i, j, k, h, g]
# Check for adjoint and bank val
if adjoint:
# Get nu-fission cross section from cell, transposed!
nfissxs = self._nfissxs[i, j, k, g, h]
else:
# Get nu-fission cross section from cell
nfissxs = self._nfissxs[i, j, k, h, g]
# set as value to be recorded
# Set as value to be recorded
val = nfissxs
# record value in matrix
@ -1312,9 +1392,11 @@ class CMFDRun(object):
i = spatial_idx[0][0]
j = spatial_idx[1][0]
k = spatial_idx[2][0]
return i, j, k, g
def _indices_to_matrix(self, i, j, k, g, ng):
# Get matrix index from coremap
matidx = ng*(self._coremap[i,j,k]) + g
return matidx
@ -1353,7 +1435,9 @@ class CMFDRun(object):
prod = sparse.csr_matrix(prod)
loss = sparse.csr_matrix(loss)
# Begin power iteration
for i in range(maxits):
# Check if reach max number of iterations
if i == maxits - 1:
raise OpenMCError('Reached maximum iterations in CMFD power '
'iteration solver.')
@ -1415,8 +1499,9 @@ class CMFDRun(object):
def _set_coremap(self):
self._mat_dim = np.sum(self._coremap)
# Create a temporary array that aggregates cumulative sum over
# accelerated regions
# Define coremap as cumulative sum over accelerated regions,
# otherwise set value to _CMFD_NOACCEL
# TODO, does algorithm work if CMFD_ACCEL is fixed number
self._coremap = np.where(self._coremap==0, _CMFD_NOACCEL,
np.cumsum(self._coremap)-1)
@ -1597,6 +1682,7 @@ class CMFDRun(object):
# Get openmc k-effective
keff = openmc.capi.keff()[0]
# Define leakage in each mesh cell and energy group
leakage = ((self._current[:,:,:,:,_CURRENTS['out_right']] - \
self._current[:,:,:,:,_CURRENTS['in_right']]) - \
(self._current[:,:,:,:,_CURRENTS['in_left']] - \
@ -1635,6 +1721,7 @@ class CMFDRun(object):
np.count_nonzero(self._resnb)))
def _compute_dtilde2(self):
#TODO add coments for this method
dtilde2 = np.zeros(self._dtilde.shape)
is_accel = self._coremap != _CMFD_NOACCEL
@ -1808,10 +1895,6 @@ class CMFDRun(object):
(self._hxyz[2] * self._diffcof[:,:,:-1,:] + \
self._hxyz[2] * neig_dc[:,:,:-1,:])), 0.0)
print("After dtilde")
print(dtilde2)
sys.exit()
def _compute_dtilde(self):
# Get maximum of spatial and group indices
nx = self._indices[0]
@ -1884,6 +1967,7 @@ class CMFDRun(object):
self._dtilde[i, j, k, g, l] = dtilde
def _compute_dhat2(self):
# TODO Write comments for this function
dhat2 = np.zeros(self._dhat.shape)
net_current_minusx = ((self._current[:,:,:,:,_CURRENTS['in_left']] - \

View file

@ -7,6 +7,8 @@ module cmfd_solver
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
use matrix_header, only: Matrix
use vector_header, only: Vector
use simulation_header, only: current_batch
use string, only: to_str
implicit none
private
@ -101,8 +103,7 @@ contains
use constants, only: ONE, ZERO
use cmfd_header, only: cmfd_shift, cmfd_ktol, cmfd_stol, cmfd_write_matrices
use simulation_header, only: keff, current_batch
use string, only: to_str
use simulation_header, only: keff
logical, intent(in) :: adjoint
@ -147,8 +148,8 @@ contains
call loss % assemble()
call prod % assemble()
if (cmfd_write_matrices) then
call loss % write('loss' // trim(to_str(current_batch)) // '.dat')
call prod % write('prod' // trim(to_str(current_batch)) // '.dat')
call loss % write('loss_gen' // trim(to_str(current_batch)) // '.dat')
call prod % write('prod_gen' // trim(to_str(current_batch)) // '.dat')
end if
! Set norms to 0
@ -176,8 +177,8 @@ contains
! Write out matrix in binary file (debugging)
if (cmfd_write_matrices) then
call loss % write('adj_loss.dat')
call prod % write('adj_prod.dat')
call loss % write('adj_loss_gen' // trim(to_str(current_batch)) // '.dat')
call prod % write('adj_prod_gen' // trim(to_str(current_batch)) // '.dat')
end if
end subroutine compute_adjoint
@ -731,17 +732,19 @@ contains
cmfd%phi = cmfd%phi/sqrt(sum(cmfd%phi*cmfd%phi))
end if
print *, cmfd % phi
! Save dominance ratio
cmfd % dom(current_batch) = norm_n/norm_o
! Write out results
if (cmfd_write_matrices) then
if (adjoint_calc) then
filename = 'adj_fluxvec.dat'
filename = 'adj_fluxvec_gen' // trim(to_str(current_batch)) // '.dat'
else
filename = 'fluxvec.dat'
filename = 'fluxvec_gen' // trim(to_str(current_batch)) // '.dat'
end if
! TODO: call phi_n % write(filename)
call phi_n % write(filename)
end if
end subroutine extract_results

View file

@ -14,7 +14,7 @@ module vector_header
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
! TODO: procedure :: write => vector_write
procedure :: write => vector_write
end type Vector
contains
@ -88,4 +88,27 @@ contains
end subroutine vector_copy
!===============================================================================
! VECTOR_WRITE write a vector to file
!===============================================================================
subroutine vector_write(self, filename)
class(Vector), target, intent(inout) :: self ! vector instance
character(*), intent(in) :: filename ! filename to output to
integer :: unit_
integer :: i
open(newunit=unit_, file=filename)
do i = 1, self % n
write(unit_,*) i, self % data(i)
print *, self % data(i)
end do
close(unit_)
end subroutine vector_write
end module vector_header