added coremap indexing

This commit is contained in:
Bryan Herman 2012-01-03 19:26:56 -08:00
parent 6f4fc51ab2
commit 51f0ab3205
2 changed files with 62 additions and 17 deletions

View file

@ -38,12 +38,11 @@ contains
! check for core map
if (allocated(cmfd % coremap)) then
call set_coremap()
print *, cmfd % coremap
end if
! compute dtilde terms
call compute_diffcoef()
STOP
! set dhats to zero
call compute_dhat()
@ -358,10 +357,10 @@ contains
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == &
& 99999 .and. cmfd % coremap(i,j,k) /= 99999) then
print *,'Before',cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)),neig_idx(1),neig_idx(2),neig_idx(3)
! get albedo
ref_albedo = get_reflector_albedo(l,g,i,j,k)
print *,'HERE --',ref_albedo
! compute dtilde
dtilde = (2*cell_dc*(1-ref_albedo))/(4*cell_dc*(1+ &
& ref_albedo)+(1-ref_albedo)*cell_hxyz(xyz_idx))
@ -805,6 +804,16 @@ use timing, only: timer_start, timer_stop
XLOOP: do i = 1,nx
! check if not including reflector
if (allocated(cmfd % coremap)) then
! check if at a reflector
if (cmfd % coremap(i,j,k) == 99999) then
cycle
end if
end if
GROUP: do g = 1,ng
! get matrix index of cell
@ -842,19 +851,30 @@ use timing, only: timer_start, timer_stop
! check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! check for core map
if (allocated(cmfd % coremap)) then
! get neighbor matrix index
neig_mat_idx = get_matrix_idx(g,neig_idx(1),neig_idx(2), &
& neig_idx(3),ng,nx,ny)
! check that neighbor is not reflector
if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= &
& 99999) then
! compute value and record to bank
val = jn/hxyz(xyz_idx)
! compute leakage coefficient for neighbor
jn = -dtilde(l) + shift_idx*dhat(l)
! record value in matrix
call MatSetValue(M,cell_mat_idx-1,neig_mat_idx-1,val, &
& INSERT_VALUES,ierr)
! get neighbor matrix index
neig_mat_idx = get_matrix_idx(g,neig_idx(1),neig_idx(2), &
& neig_idx(3),ng,nx,ny)
! compute value and record to bank
val = jn/hxyz(xyz_idx)
! record value in matrix
call MatSetValue(M,cell_mat_idx-1,neig_mat_idx-1,val, &
& INSERT_VALUES,ierr)
end if
end if
end if
@ -962,6 +982,16 @@ use timing, only: timer_start, timer_stop
XLOOP: do i = 1,nx
! check if not including reflector
if (allocated(cmfd % coremap)) then
! check if at a reflector
if (cmfd % coremap(i,j,k) == 99999) then
cycle
end if
end if
GROUP: do g = 1,ng
NFISS: do h = 1,ng

View file

@ -83,8 +83,18 @@ contains
! local variables
integer :: nidx ! index in matrix
! compute index
nidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
! check if coremap is used
if (allocated(cmfd % coremap)) then
! get idx from core map
nidx = ng*(cmfd % coremap(i,j,k)) + (ng - g)
else
! compute index
nidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1)
end if
! record value to function
get_matrix_idx = nidx
@ -428,12 +438,17 @@ contains
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
integer :: nmat ! number of fuel assy, dim of mat
! extract spatial indices from object
nx = cmfd % indices(1)
ny = cmfd % indices(2)
nz = cmfd % indices(3)
! count how many fuel assemblies exist
nmat = sum(cmfd % coremap - 1)
print *,'Dim is ',nmat
! begin loops over spatial indices
ZLOOP: do k = 1,nz
@ -490,7 +505,7 @@ contains
! calculate albedo
albedo = (current(2*l-1)/current(2*l))**(shift_idx)
print *,current(2*l-1),current(2*l)
! assign to function variable
get_reflector_albedo = albedo