mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
added module to handle the interface with openmc to compute xs and diffusion parameters
This commit is contained in:
Bryan Herman
parent
c886be7433
commit
daeb8f376c
1 changed files with 826 additions and 0 deletions
826
src/cmfd_data.F90
Normal file
826
src/cmfd_data.F90
Normal file
|
|
@ -0,0 +1,826 @@
|
|||
module cmfd_data
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!==============================================================================
|
||||
! SET_UP_CMFD
|
||||
!===============================================================================
|
||||
|
||||
subroutine set_up_cmfd()
|
||||
|
||||
use global, only: cmfd
|
||||
use cmfd_header, only: allocate_cmfd
|
||||
|
||||
! initialize data
|
||||
call allocate_cmfd(cmfd)
|
||||
|
||||
! calculate all cross sections based on reaction rates from last batch
|
||||
call compute_xs()
|
||||
|
||||
! write out the neutron balance file
|
||||
! call neutron_balance()
|
||||
|
||||
! check for core map
|
||||
if (allocated(cmfd % coremap)) then
|
||||
call set_coremap()
|
||||
end if
|
||||
|
||||
! compute dtilde terms
|
||||
call compute_diffcoef()
|
||||
|
||||
! set dhats to zero
|
||||
call compute_dhat()
|
||||
|
||||
end subroutine set_up_cmfd
|
||||
|
||||
!===============================================================================
|
||||
! READ_INPUT reads the CMFD input file and organizes it into a data structure
|
||||
!===============================================================================
|
||||
|
||||
subroutine read_cmfd_xml()
|
||||
|
||||
use global
|
||||
use string
|
||||
use xml_data_cmfd_t
|
||||
|
||||
integer :: ng=1 ! number of energy groups (default 1)
|
||||
integer :: n_words ! number of words read
|
||||
logical :: file_exists ! does cmfd.xml exist?
|
||||
character(MAX_LINE_LEN) :: filename
|
||||
character(MAX_WORD_LEN) :: words(MAX_WORDS)
|
||||
|
||||
! read cmfd infput file
|
||||
filename = "cmfd.xml"
|
||||
inquire(FILE=filename, EXIST=file_exists)
|
||||
if (.not. file_exists) then
|
||||
write(*,*) "Cannot perform CMFD"
|
||||
STOP
|
||||
end if
|
||||
|
||||
! parse cmfd.xml file
|
||||
call read_xml_file_cmfd_t(filename)
|
||||
|
||||
! set spatial dimensions in cmfd object
|
||||
cmfd % indices(1:3) = mesh_ % dimension(1:3) ! sets spatial dimensions
|
||||
|
||||
! get number of energy groups
|
||||
if (len_trim(mesh_ % energy) > 0) then
|
||||
call split_string(mesh_ % energy, words, n_words)
|
||||
ng = n_words - 1
|
||||
end if
|
||||
cmfd % indices(4) = ng ! sets energy group dimension
|
||||
|
||||
! set global albedo
|
||||
cmfd % albedo = mesh_ % albedo
|
||||
|
||||
! get acceleration map
|
||||
if (associated(mesh_ % map)) then
|
||||
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
|
||||
& cmfd % indices(3)))
|
||||
cmfd % coremap = reshape(mesh_ % map,(cmfd % indices(1:3)))
|
||||
end if
|
||||
|
||||
! check for core map activation by printing note
|
||||
if (allocated(cmfd % coremap)) print *,"Core Map Overlay Activated"
|
||||
|
||||
! create tally objects
|
||||
call create_cmfd_tally()
|
||||
|
||||
end subroutine read_cmfd_xml
|
||||
|
||||
!===============================================================================
|
||||
! CREATE_CMFD_TALLY creates the tally object for OpenMC to process for CMFD
|
||||
! accleration.
|
||||
! There are 3 tally types:
|
||||
! 1: Only an energy in filter-> flux,total,p1 scatter
|
||||
! 2: Energy in and energy out filter-> nu-scatter,nu-fission
|
||||
! 3: Surface current
|
||||
!===============================================================================
|
||||
|
||||
subroutine create_cmfd_tally()
|
||||
|
||||
use datatypes, only: dict_add_key, dict_get_key
|
||||
use error, only: fatal_error, warning
|
||||
use global
|
||||
use mesh_header, only: StructuredMesh
|
||||
use string
|
||||
use tally_header, only: TallyObject, TallyScore
|
||||
use xml_data_cmfd_t
|
||||
|
||||
integer :: i ! loop counter
|
||||
integer :: j ! loop counter
|
||||
integer :: id ! user-specified identifier
|
||||
integer :: index ! index in mesh array
|
||||
integer :: n ! size of arrays in mesh specification
|
||||
integer :: ng=1 ! number of energy groups (default 1)
|
||||
integer :: n_words ! number of words read
|
||||
character(MAX_LINE_LEN) :: filename
|
||||
character(MAX_WORD_LEN) :: words(MAX_WORDS)
|
||||
type(TallyObject), pointer :: t => null()
|
||||
type(StructuredMesh), pointer :: m => null()
|
||||
|
||||
! parse cmfd.xml file
|
||||
filename = trim(path_input) // "cmfd.xml"
|
||||
call read_xml_file_cmfd_t(filename)
|
||||
|
||||
! allocate mesh
|
||||
n_meshes = 1
|
||||
allocate(meshes(n_meshes))
|
||||
m => meshes(1)
|
||||
|
||||
! set mesh id
|
||||
m % id = 1
|
||||
|
||||
! set mesh type to rectangular
|
||||
m % type = LATTICE_RECT
|
||||
|
||||
! determine number of dimensions for mesh
|
||||
n = size(mesh_ % dimension)
|
||||
if (n /= 2 .and. n /= 3) then
|
||||
message = "Mesh must be two or three dimensions."
|
||||
call fatal_error()
|
||||
end if
|
||||
m % n_dimension = n
|
||||
|
||||
! allocate attribute arrays
|
||||
allocate(m % dimension(n))
|
||||
allocate(m % origin(n))
|
||||
allocate(m % width(n))
|
||||
allocate(m % upper_right(n))
|
||||
|
||||
! read dimensions in each direction
|
||||
m % dimension = mesh_ % dimension
|
||||
|
||||
! read mesh origin location
|
||||
if (m % n_dimension /= size(mesh_ % origin)) then
|
||||
message = "Number of entries on <origin> must be the same as " // &
|
||||
"the number of entries on <dimension>."
|
||||
call fatal_error()
|
||||
end if
|
||||
m % origin = mesh_ % origin
|
||||
|
||||
! read mesh widths
|
||||
if (size(mesh_ % width) /= size(mesh_ % origin)) then
|
||||
message = "Number of entries on <width> must be the same as " // &
|
||||
"the number of entries on <origin>."
|
||||
call fatal_error()
|
||||
end if
|
||||
m % width = mesh_ % width
|
||||
|
||||
! set upper right coordinate
|
||||
m % upper_right = m % origin + m % dimension * m % width
|
||||
|
||||
! add mesh to dictionary
|
||||
call dict_add_key(mesh_dict, m % id, 1)
|
||||
|
||||
! allocate tallies
|
||||
n_tallies = 3
|
||||
allocate(tallies(n_tallies))
|
||||
|
||||
! begin loop around tallies
|
||||
do i = 1,n_tallies
|
||||
t => tallies(i)
|
||||
|
||||
! allocate arrays for number of bins and stride in scores array
|
||||
allocate(t % n_bins(TALLY_TYPES))
|
||||
allocate(t % stride(TALLY_TYPES))
|
||||
|
||||
! initialize number of bins and stride
|
||||
t % n_bins = 0
|
||||
t % stride = 0
|
||||
|
||||
! record tally id which is equivalent to loop number
|
||||
t % id = i
|
||||
|
||||
! set mesh filter mesh id = 1
|
||||
t % mesh = 1
|
||||
m => meshes(1)
|
||||
t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
|
||||
|
||||
! read and set incoming energy mesh filter
|
||||
if (len_trim(mesh_ % energy) > 0) then
|
||||
call split_string(mesh_ % energy,words,n_words)
|
||||
ng = n_words
|
||||
allocate(t % energy_in(n_words))
|
||||
do j = 1,n_words
|
||||
t % energy_in(j) = str_to_real(words(j))
|
||||
end do
|
||||
t % n_bins(T_ENERGYIN) = n_words - 1
|
||||
end if
|
||||
|
||||
if (i == 1) then
|
||||
|
||||
! allocate macro reactions
|
||||
allocate(t % macro_bins(3))
|
||||
t % n_macro_bins = 3
|
||||
|
||||
! set macro_bins
|
||||
t % macro_bins(1) % scalar = MACRO_FLUX
|
||||
t % macro_bins(2) % scalar = MACRO_TOTAL
|
||||
t % macro_bins(3) % scalar = MACRO_SCATTER_1
|
||||
|
||||
else if (i == 2) then
|
||||
|
||||
! read and set outgoing energy mesh filter
|
||||
if (len_trim(mesh_ % energy) > 0) then
|
||||
call split_string(mesh_ % energy, words, n_words)
|
||||
allocate(t % energy_out(n_words))
|
||||
do j = 1, n_words
|
||||
t % energy_out(j) = str_to_real(words(j))
|
||||
end do
|
||||
t % n_bins(T_ENERGYOUT) = n_words - 1
|
||||
end if
|
||||
|
||||
! allocate macro reactions
|
||||
allocate(t % macro_bins(2))
|
||||
t % n_macro_bins = 2
|
||||
|
||||
! set macro_bins
|
||||
t % macro_bins(1) % scalar = MACRO_NU_SCATTER
|
||||
t % macro_bins(2) % scalar = MACRO_NU_FISSION
|
||||
|
||||
else if (i == 3) then
|
||||
|
||||
! allocate macro reactions
|
||||
allocate(t % macro_bins(1))
|
||||
t % n_macro_bins = 1
|
||||
|
||||
! set macro bins
|
||||
t % macro_bins(1) % scalar = MACRO_CURRENT
|
||||
t % surface_current = .true.
|
||||
|
||||
! since the number of bins for the mesh filter was already set
|
||||
! assuming it was a flux tally, we need to adjust the number of
|
||||
! bins
|
||||
t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
|
||||
|
||||
! get pointer to mesh
|
||||
id = t % mesh
|
||||
index = dict_get_key(mesh_dict, id)
|
||||
m => meshes(index)
|
||||
|
||||
! we need to increase the dimension by one since we also need
|
||||
! currents coming into and out of the boundary mesh cells.
|
||||
if (size(m % dimension) == 2) then
|
||||
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
|
||||
& product(m % dimension + 1) * 4
|
||||
elseif (size(m % dimension) == 3) then
|
||||
t % n_bins(T_MESH) = t % n_bins(T_MESH) + &
|
||||
product(m % dimension + 1) * 6
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
end subroutine create_cmfd_tally
|
||||
|
||||
!===============================================================================
|
||||
! COMPUTE_XS takes tallies and computes macroscopic cross sections
|
||||
!===============================================================================
|
||||
|
||||
subroutine compute_xs()
|
||||
|
||||
use global
|
||||
use mesh, only: mesh_indices_to_bin
|
||||
use mesh_header, only: StructuredMesh
|
||||
use tally_header, only: TallyObject, TallyScore
|
||||
|
||||
integer :: nx ! number of mesh cells in x direction
|
||||
integer :: ny ! number of mesh cells in y direction
|
||||
integer :: nz ! number of mesh cells in z direction
|
||||
integer :: ng ! number of energy groups
|
||||
integer :: i ! iteration counter for x
|
||||
integer :: j ! iteration counter for y
|
||||
integer :: k ! iteration counter for z
|
||||
integer :: g ! iteration counter for g
|
||||
integer :: h ! iteration counter for outgoing groups
|
||||
integer :: ijk(3) ! indices for mesh cell
|
||||
integer :: score_index ! index to pull from tally object
|
||||
integer :: bins(TALLY_TYPES) ! bins for filters
|
||||
|
||||
real(8) :: flux ! temp variable for flux
|
||||
|
||||
type(TallyObject), pointer :: t ! pointer for tally object
|
||||
type(StructuredMesh), pointer :: m ! pointer for mesh object
|
||||
|
||||
! extract spatial and energy indices from object
|
||||
nx = cmfd % indices(1)
|
||||
ny = cmfd % indices(2)
|
||||
nz = cmfd % indices(3)
|
||||
ng = cmfd % indices(4)
|
||||
|
||||
! begin loop around space and energy groups
|
||||
ZLOOP: do k = 1,nz
|
||||
|
||||
YLOOP: do j = 1,ny
|
||||
|
||||
XLOOP: do i = 1,nx
|
||||
|
||||
OUTGROUP: do h = 1,ng
|
||||
|
||||
! begin with first tally
|
||||
t => tallies(1)
|
||||
m => meshes(t % mesh)
|
||||
|
||||
! set mesh widths
|
||||
cmfd % hxyz(1,:,:,:) = m % width(1) ! set x width
|
||||
cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width
|
||||
cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width
|
||||
|
||||
! reset all bins to 1
|
||||
bins = 1
|
||||
|
||||
! set ijk as mesh indices
|
||||
ijk = (/ i, j, k /)
|
||||
|
||||
! get bin number for mesh indices
|
||||
bins(T_MESH) = mesh_indices_to_bin(m,ijk)
|
||||
|
||||
! apply energy in filter
|
||||
bins(T_ENERGYIN) = ng - h + 1
|
||||
|
||||
! calculate score index from bins
|
||||
score_index = sum((bins - 1) * t%stride) + 1
|
||||
|
||||
! get flux
|
||||
flux = t % scores(score_index,1) % val
|
||||
cmfd % flux(h,i,j,k) = flux
|
||||
|
||||
! detect zero flux
|
||||
if ((flux - 0.0D0) < 1.0D-10) then
|
||||
if (.not. allocated(cmfd%coremap)) then
|
||||
write(*,*) 'Fatal: detected zero flux without coremap'
|
||||
stop
|
||||
else
|
||||
write(*,*) 'Warning: detected zero flux at:',i,j,k
|
||||
flux = 99999.0D0
|
||||
if (.not. cmfd%coremap(i,j,k) == 99999) then
|
||||
write(*,*) 'Fatal: need to check core map with zero flux'
|
||||
stop
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
|
||||
! get total rr and convert to total xs
|
||||
cmfd % totalxs(h,i,j,k) = t % scores(score_index,2) % val / flux
|
||||
|
||||
! get p1 scatter rr and convert to p1 scatter xs
|
||||
cmfd % p1scattxs(h,i,j,k) = t % scores(score_index,3) % val / flux
|
||||
|
||||
! calculate diffusion coefficient
|
||||
cmfd % diffcof(h,i,j,k) = 1/(3*(cmfd % totalxs(h,i,j,k) - &
|
||||
& cmfd % p1scattxs(h,i,j,k)))
|
||||
cmfd % diffcof(h,i,j,k) = 1/(3*(cmfd % totalxs(h,i,j,k)))
|
||||
|
||||
! begin loop to get energy out tallies
|
||||
INGROUP: do g = 1,ng
|
||||
|
||||
! associate tally pointer to energy out tally object
|
||||
t => tallies(2)
|
||||
|
||||
! set energy out bin
|
||||
bins(T_ENERGYOUT) = ng - g + 1
|
||||
|
||||
! calculate score index from bins
|
||||
score_index = sum((bins - 1) * t%stride) + 1
|
||||
|
||||
! get scattering
|
||||
cmfd % scattxs(h,g,i,j,k) = t % scores(score_index,1) % val / flux
|
||||
|
||||
! get nu-fission
|
||||
cmfd % nfissxs(h,g,i,j,k) = t % scores(score_index,2) % val / flux
|
||||
|
||||
end do INGROUP
|
||||
|
||||
! extract surface currents
|
||||
t => tallies(3)
|
||||
|
||||
! initialize and filter for energy
|
||||
bins = 1
|
||||
bins(TS_ENERGYIN) = ng - h + 1
|
||||
|
||||
! left surface
|
||||
bins(1:3) = (/ i-1, j, k /) + 1
|
||||
bins(TS_SURFACE) = IN_RIGHT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(1,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_RIGHT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(2,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
! right surface
|
||||
bins(1:3) = (/ i, j, k /) + 1
|
||||
bins(TS_SURFACE) = IN_RIGHT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(3,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_RIGHT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(4,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
! back surface
|
||||
bins(1:3) = (/ i, j-1, k /) + 1
|
||||
bins(TS_SURFACE) = IN_FRONT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(5,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_FRONT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(6,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
! front surface
|
||||
bins(1:3) = (/ i, j, k /) + 1
|
||||
bins(TS_SURFACE) = IN_FRONT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(7,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_FRONT
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(8,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
! bottom surface
|
||||
bins(1:3) = (/ i, j, k-1 /) + 1
|
||||
bins(TS_SURFACE) = IN_TOP
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(9,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_TOP
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(10,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
! top surface
|
||||
bins(1:3) = (/ i, j, k /) + 1
|
||||
bins(TS_SURFACE) = IN_TOP
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! incoming
|
||||
cmfd % current(11,h,i,j,k) = t % scores(score_index,1) % val
|
||||
bins(TS_SURFACE) = OUT_TOP
|
||||
score_index = sum((bins - 1) * t % stride) + 1 ! outgoing
|
||||
cmfd % current(12,h,i,j,k) = t % scores(score_index,1) % val
|
||||
|
||||
end do OUTGROUP
|
||||
|
||||
end do XLOOP
|
||||
|
||||
end do YLOOP
|
||||
|
||||
end do ZLOOP
|
||||
|
||||
end subroutine compute_xs
|
||||
|
||||
!===============================================================================
|
||||
! COMPUTE_DIFFCOEF computes the diffusion coupling coefficient
|
||||
!===============================================================================
|
||||
|
||||
subroutine compute_diffcoef()
|
||||
|
||||
use global, only: cmfd
|
||||
|
||||
! local variables
|
||||
integer :: nx ! maximum number of cells in x direction
|
||||
integer :: ny ! maximum number of cells in y direction
|
||||
integer :: nz ! maximum number of cells in z direction
|
||||
integer :: ng ! maximum number of energy groups
|
||||
integer :: nxyz(3,2) ! single vector containing boundary locations
|
||||
integer :: i ! iteration counter for x
|
||||
integer :: j ! iteration counter for y
|
||||
integer :: k ! iteration counter for z
|
||||
integer :: g ! iteration counter for groups
|
||||
integer :: l ! iteration counter for leakages
|
||||
integer :: xyz_idx ! index for determining if x,y or z leakage
|
||||
integer :: dir_idx ! index for determining - or + face of cell
|
||||
integer :: shift_idx ! parameter to shift index by +1 or -1
|
||||
integer :: neig_idx(3) ! spatial indices of neighbour
|
||||
integer :: bound(6) ! vector containing indices for boudary check
|
||||
real(8) :: albedo(6) ! albedo vector with global boundaries
|
||||
real(8) :: cell_totxs ! total cross section of current ijk cell
|
||||
real(8) :: cell_dc ! diffusion coef of current cell
|
||||
real(8) :: cell_hxyz(3) ! cell dimensions of current ijk cell
|
||||
real(8) :: neig_totxs ! total xs of neighbor cell
|
||||
real(8) :: neig_dc ! diffusion coefficient of neighbor cell
|
||||
real(8) :: neig_hxyz(3) ! cell dimensions of neighbor cell
|
||||
real(8) :: dtilde ! finite difference coupling parameter
|
||||
real(8) :: ref_albedo ! albedo to reflector
|
||||
|
||||
! get maximum of spatial and group indices
|
||||
nx = cmfd%indices(1)
|
||||
ny = cmfd%indices(2)
|
||||
nz = cmfd%indices(3)
|
||||
ng = cmfd%indices(4)
|
||||
|
||||
! create single vector of these indices for boundary calculation
|
||||
nxyz(1,:) = (/1,nx/)
|
||||
nxyz(2,:) = (/1,ny/)
|
||||
nxyz(3,:) = (/1,nz/)
|
||||
|
||||
! get boundary condition information
|
||||
albedo = cmfd%albedo
|
||||
|
||||
! geting loop over group and spatial indices
|
||||
ZLOOP: do k = 1,nz
|
||||
|
||||
YLOOP: do j = 1,ny
|
||||
|
||||
XLOOP: do i = 1,nx
|
||||
|
||||
GROUP: do g = 1,ng
|
||||
|
||||
! get cell data
|
||||
cell_dc = cmfd%diffcof(g,i,j,k)
|
||||
cell_hxyz = cmfd%hxyz(:,i,j,k)
|
||||
|
||||
|
||||
! setup of vector to identify boundary conditions
|
||||
bound = (/i,i,j,j,k,k/)
|
||||
|
||||
! begin loop around sides of cell for leakage
|
||||
LEAK: do l = 1,6
|
||||
|
||||
! define xyz and +/- indices
|
||||
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
|
||||
dir_idx = 2 - mod(l,2) ! -=1, +=2
|
||||
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
|
||||
|
||||
! check if at a boundary
|
||||
if (bound(l) == nxyz(xyz_idx,dir_idx)) then
|
||||
|
||||
! compute dtilde
|
||||
dtilde = (2*cell_dc*(1-albedo(l)))/(4*cell_dc*(1+albedo(l)) + &
|
||||
& (1-albedo(l))*cell_hxyz(xyz_idx))
|
||||
|
||||
else ! not a boundary
|
||||
|
||||
! compute neighboring cell indices
|
||||
neig_idx = (/i,j,k/) ! begin with i,j,k
|
||||
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
|
||||
|
||||
! get neigbor cell data
|
||||
neig_dc = cmfd%diffcof(g,neig_idx(1),neig_idx(2),neig_idx(3))
|
||||
neig_hxyz = cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3))
|
||||
|
||||
! check for fuel-reflector interface
|
||||
if (allocated(cmfd % coremap)) then
|
||||
|
||||
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == &
|
||||
& 99999 .and. cmfd % coremap(i,j,k) /= 99999) then
|
||||
|
||||
! get albedo
|
||||
ref_albedo = get_reflector_albedo(l,g,i,j,k)
|
||||
|
||||
! compute dtilde
|
||||
dtilde = (2*cell_dc*(1-ref_albedo))/(4*cell_dc*(1+ &
|
||||
& ref_albedo)+(1-ref_albedo)*cell_hxyz(xyz_idx))
|
||||
|
||||
else ! not next to a reflector or no core map
|
||||
|
||||
! compute dtilde
|
||||
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
|
||||
& cell_hxyz(xyz_idx)*neig_dc)
|
||||
|
||||
end if
|
||||
|
||||
else ! no core map
|
||||
|
||||
! compute dtilde
|
||||
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
|
||||
& cell_hxyz(xyz_idx)*neig_dc)
|
||||
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
! record dtilde in cmfd object
|
||||
cmfd%dtilde(l,g,i,j,k) = dtilde
|
||||
|
||||
end do LEAK
|
||||
|
||||
end do GROUP
|
||||
|
||||
end do XLOOP
|
||||
|
||||
end do YLOOP
|
||||
|
||||
end do ZLOOP
|
||||
|
||||
end subroutine compute_diffcoef
|
||||
|
||||
!===============================================================================
|
||||
! COMPUTE_DHAT computes the nonlinear coupling coefficient
|
||||
!===============================================================================
|
||||
|
||||
subroutine compute_dhat()
|
||||
|
||||
use global, only:cmfd
|
||||
|
||||
! local variables
|
||||
integer :: nx ! maximum number of cells in x direction
|
||||
integer :: ny ! maximum number of cells in y direction
|
||||
integer :: nz ! maximum number of cells in z direction
|
||||
integer :: ng ! maximum number of energy groups
|
||||
integer :: nxyz(3,2) ! single vector containing boundary locations
|
||||
integer :: i ! iteration counter for x
|
||||
integer :: j ! iteration counter for y
|
||||
integer :: k ! iteration counter for z
|
||||
integer :: g ! iteration counter for groups
|
||||
integer :: l ! iteration counter for leakages
|
||||
integer :: xyz_idx ! index for determining if x,y or z leakage
|
||||
integer :: dir_idx ! index for determining - or + face of cell
|
||||
integer :: shift_idx ! parameter to shift index by +1 or -1
|
||||
integer :: neig_idx(3) ! spatial indices of neighbour
|
||||
integer :: bound(6) ! vector containing indices for boudary check
|
||||
real(8) :: cell_dtilde(6) ! cell dtilde for each face
|
||||
real(8) :: cell_flux ! flux in current cell
|
||||
real(8) :: current(12) ! area integrated cell current at each face
|
||||
real(8) :: net_current ! net current on a face
|
||||
real(8) :: neig_flux ! flux in neighbor cell
|
||||
real(8) :: dhat ! dhat equivalence parameter
|
||||
|
||||
! get maximum of spatial and group indices
|
||||
nx = cmfd%indices(1)
|
||||
ny = cmfd%indices(2)
|
||||
nz = cmfd%indices(3)
|
||||
ng = cmfd%indices(4)
|
||||
|
||||
! create single vector of these indices for boundary calculation
|
||||
nxyz(1,:) = (/1,nx/)
|
||||
nxyz(2,:) = (/1,ny/)
|
||||
nxyz(3,:) = (/1,nz/)
|
||||
|
||||
! geting loop over group and spatial indices
|
||||
ZLOOP: do k = 1,nz
|
||||
|
||||
YLOOP: do j = 1,ny
|
||||
|
||||
XLOOP: do i = 1,nx
|
||||
|
||||
GROUP: do g = 1,ng
|
||||
|
||||
! get cell data
|
||||
cell_dtilde = cmfd%dtilde(:,g,i,j,k)
|
||||
cell_flux = cmfd%flux(g,i,j,k)/product(cmfd%hxyz(:,i,j,k))
|
||||
current = cmfd%current(:,g,i,j,k)
|
||||
|
||||
! setup of vector to identify boundary conditions
|
||||
bound = (/i,i,j,j,k,k/)
|
||||
|
||||
! begin loop around sides of cell for leakage
|
||||
LEAK: do l = 1,6
|
||||
|
||||
! define xyz and +/- indices
|
||||
xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3
|
||||
dir_idx = 2 - mod(l,2) ! -=1, +=2
|
||||
shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1
|
||||
|
||||
! calculate net current on l face (divided by surf area)
|
||||
net_current = (current(2*l) - current(2*l-1)) / &
|
||||
& product(cmfd%hxyz(:,i,j,k)) * &
|
||||
& cmfd%hxyz(xyz_idx,i,j,k)
|
||||
|
||||
! check if at a boundary
|
||||
if (bound(l) == nxyz(xyz_idx,dir_idx)) then
|
||||
|
||||
! compute dhat
|
||||
dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) / &
|
||||
& cell_flux
|
||||
|
||||
else ! not a boundary
|
||||
|
||||
! compute neighboring cell indices
|
||||
neig_idx = (/i,j,k/) ! begin with i,j,k
|
||||
neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx)
|
||||
|
||||
! get neigbor flux
|
||||
neig_flux = cmfd%flux(g,neig_idx(1),neig_idx(2),neig_idx(3)) / &
|
||||
product(cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3)))
|
||||
|
||||
! check for fuel-reflector interface
|
||||
if (allocated(cmfd % coremap)) then
|
||||
|
||||
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == &
|
||||
& 99999 .and. cmfd % coremap(i,j,k) /= 99999) then
|
||||
|
||||
! compute dhat
|
||||
dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) /&
|
||||
& cell_flux
|
||||
|
||||
else ! not a fuel-reflector interface
|
||||
|
||||
! compute dhat
|
||||
dhat = (net_current + shift_idx*cell_dtilde(l)* &
|
||||
& (neig_flux - cell_flux))/(neig_flux + cell_flux)
|
||||
|
||||
end if
|
||||
|
||||
else ! not for fuel-reflector case
|
||||
|
||||
! compute dhat
|
||||
dhat = (net_current + shift_idx*cell_dtilde(l)* &
|
||||
& (neig_flux - cell_flux))/(neig_flux + cell_flux)
|
||||
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
! record dtilde in cmfd object
|
||||
cmfd%dhat(l,g,i,j,k) = dhat
|
||||
|
||||
end do LEAK
|
||||
|
||||
end do GROUP
|
||||
|
||||
end do XLOOP
|
||||
|
||||
end do YLOOP
|
||||
|
||||
end do ZLOOP
|
||||
|
||||
end subroutine compute_dhat
|
||||
|
||||
!===============================================================================
|
||||
! SET_COREMAP is a routine that sets the core mapping information
|
||||
!===============================================================================
|
||||
|
||||
subroutine set_coremap()
|
||||
|
||||
use global, only: cmfd
|
||||
|
||||
integer :: kount=1 ! counter for unique fuel assemblies
|
||||
integer :: nx ! number of mesh cells in x direction
|
||||
integer :: ny ! number of mesh cells in y direction
|
||||
integer :: nz ! number of mesh cells in z direction
|
||||
integer :: ng ! number of energy groups
|
||||
integer :: i ! iteration counter for x
|
||||
integer :: j ! iteration counter for y
|
||||
integer :: k ! iteration counter for z
|
||||
|
||||
! extract spatial indices from object
|
||||
nx = cmfd % indices(1)
|
||||
ny = cmfd % indices(2)
|
||||
nz = cmfd % indices(3)
|
||||
|
||||
! count how many fuel assemblies exist
|
||||
cmfd % mat_dim = sum(cmfd % coremap - 1)
|
||||
|
||||
! begin loops over spatial indices
|
||||
ZLOOP: do k = 1,nz
|
||||
|
||||
YLOOP: do j = 1,ny
|
||||
|
||||
XLOOP: do i = 1,nx
|
||||
|
||||
! check for reflector
|
||||
if (cmfd % coremap(i,j,k) == 1) then
|
||||
|
||||
! reset value to 99999
|
||||
cmfd % coremap(i,j,k) = 99999
|
||||
|
||||
else
|
||||
|
||||
! must be a fuel --> give unique id number
|
||||
cmfd % coremap(i,j,k) = kount
|
||||
kount = kount + 1
|
||||
|
||||
end if
|
||||
|
||||
end do XLOOP
|
||||
|
||||
end do YLOOP
|
||||
|
||||
end do ZLOOP
|
||||
|
||||
end subroutine set_coremap
|
||||
|
||||
!===============================================================================
|
||||
! GET_REFLECTOR_ALBEDO is a function that calculates the albedo to the reflector
|
||||
!===============================================================================
|
||||
|
||||
function get_reflector_albedo(l,g,i,j,k)
|
||||
|
||||
use global, only: cmfd
|
||||
|
||||
! function variable
|
||||
real(8) :: get_reflector_albedo ! reflector albedo
|
||||
|
||||
! local variable
|
||||
integer :: i ! iteration counter for x
|
||||
integer :: j ! iteration counter for y
|
||||
integer :: k ! iteration counter for z
|
||||
integer :: g ! iteration counter for groups
|
||||
integer :: l ! iteration counter for leakages
|
||||
integer :: shift_idx ! parameter to shift index by +1 or -1
|
||||
real(8) :: current(12) ! partial currents for all faces of mesh cell
|
||||
real(8) :: albedo ! the albedo
|
||||
|
||||
! get partial currents from object
|
||||
current = cmfd%current(:,g,i,j,k)
|
||||
|
||||
! define xyz and +/- indices
|
||||
shift_idx = -2*mod(l,2) + 1 ! shift neig by -1 or +1
|
||||
|
||||
! calculate albedo
|
||||
albedo = (current(2*l-1)/current(2*l))**(shift_idx)
|
||||
|
||||
! assign to function variable
|
||||
get_reflector_albedo = albedo
|
||||
|
||||
end function get_reflector_albedo
|
||||
|
||||
end module cmfd_data
|
||||
Loading…
Add table
Add a link
Reference in a new issue