diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 3b46e6de69..cb92f708e9 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -21,11 +21,18 @@ contains use cmfd_header, only: allocate_cmfd use constants, only: CMFD_NOACCEL use global, only: cmfd, cmfd_coremap, cmfd_downscatter, & - cmfd_fix_balance + cmfd_fix_balance, message + use output, only: write_message + + message = 'Extracting CMFD data...' + call write_message(1) ! Check for core map and set it up if ((cmfd_coremap) .and. (cmfd%mat_dim == CMFD_NOACCEL)) call set_coremap() + ! Extract tallies to put in cmfd data object + call extract_tallies() + ! Calculate all cross sections based on reaction rates from last batch call compute_xs() @@ -47,10 +54,10 @@ contains end subroutine set_up_cmfd !=============================================================================== -! COMPUTE_XS takes tallies and computes macroscopic cross sections +! EXTRACT_TALLIES !=============================================================================== - subroutine compute_xs() + subroutine extract_tallies() use constants, only: FILTER_MESH, FILTER_ENERGYIN, FILTER_ENERGYOUT, & FILTER_SURFACE, IN_RIGHT, OUT_RIGHT, IN_FRONT, & @@ -58,7 +65,7 @@ contains ONE, TINY_BIT use error, only: fatal_error use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,& - matching_bins, current_batch + matching_bins, current_batch, cmfd_n_save use mesh, only: mesh_indices_to_bin use mesh_header, only: StructuredMesh use string, only: to_str @@ -73,7 +80,7 @@ contains integer :: k ! iteration counter for z integer :: g ! iteration counter for g integer :: h ! iteration counter for outgoing groups - integer :: l ! iteration counter for leakage debug + integer :: b ! current batch index in moving window integer :: ital ! tally object index integer :: ijk(3) ! indices for mesh cell integer :: score_index ! index to pull from tally object @@ -82,12 +89,9 @@ contains integer :: i_filter_ein ! index for incoming energy filter integer :: i_filter_eout ! index for outgoing energy filter integer :: i_filter_surf ! index for surface filter - real(8) :: flux ! temp variable for flux - real(8) :: leak1 ! group 1 leakage - real(8) :: leak2 ! group 2 leakage type(TallyObject), pointer :: t => null() ! pointer for tally object type(StructuredMesh), pointer :: m => null() ! pointer for mesh object - +print *, 'CMFD INDEX:', cmfd % idx ! Extract spatial and energy indices from object nx = cmfd % indices(1) ny = cmfd % indices(2) @@ -108,6 +112,10 @@ contains cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width + ! Set batch index for moving window + b = cmfd % idx + + ! Set balance keff to zero cmfd % keff_bal = ZERO ! Begin loop around tallies @@ -161,11 +169,10 @@ contains score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1 ! Get flux - flux = t % results(1,score_index) % sum - cmfd % flux(h,i,j,k) = flux + cmfd % flux_rate(h,i,j,k,b) = t % results(1,score_index) % sum ! Detect zero flux, abort if located - if ((flux - ZERO) < TINY_BIT) then + if ((cmfd % flux_rate(h,i,j,k,b) - ZERO) < TINY_BIT) then message = 'Detected zero flux without coremap overlay at: (' & // to_str(i) // ',' // to_str(j) // ',' // to_str(k) & // ') in group ' // to_str(h) @@ -173,14 +180,10 @@ contains end if ! Get total rr and convert to total xs - cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux + cmfd % total_rate(h,i,j,k,b) = t % results(2,score_index) % sum ! Get p1 scatter rr and convert to p1 scatter xs - cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux - - ! Calculate diffusion coefficient - cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - & - cmfd % p1scattxs(h,i,j,k))) + cmfd % p1scatt_rate(h,i,j,k,b) = t % results(3,score_index) % sum else if (ital == 2) then @@ -208,12 +211,10 @@ contains score_index = sum((matching_bins(1:t%n_filters) - 1) * t%stride) + 1 ! Get scattering - cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /& - cmfd % flux(h,i,j,k) + cmfd % scatt_rate(h,g,i,j,k,b) = t % results(1,score_index) % sum ! Get nu-fission - cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /& - cmfd % flux(h,i,j,k) + cmfd % nfiss_rate(h,g,i,j,k,b) = t % results(2,score_index) % sum ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & @@ -235,60 +236,60 @@ contains (/ i-1, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(1,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(2,h,i,j,k,b) = t % results(1,score_index) % sum ! Right surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(3,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(4,h,i,j,k,b) = t % results(1,score_index) % sum ! Back surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j-1, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(5,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(6,h,i,j,k,b) = t % results(1,score_index) % sum ! Front surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(7,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(8,h,i,j,k,b) = t % results(1,score_index) % sum ! Bottom surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k-1 /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(9,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(10,h,i,j,k,b) = t % results(1,score_index) % sum ! Top surface matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! incoming - cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(11,h,i,j,k,b) = t % results(1,score_index) % sum matching_bins(i_filter_surf) = OUT_TOP score_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 ! outgoing - cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum + cmfd % current_rate(12,h,i,j,k,b) = t % results(1,score_index) % sum end if TALLY @@ -309,6 +310,109 @@ contains if (associated(t)) nullify(t) if (associated(m)) nullify(m) + ! Move batch index + cmfd % idx = cmfd % idx + 1 + if (cmfd % idx > cmfd_n_save) then + cmfd % idx = 1 + end if + + end subroutine extract_tallies + +!=============================================================================== +! COMPUTE_XS takes tallies and computes macroscopic cross sections +!=============================================================================== + + subroutine compute_xs() + + use constants, only: CMFD_NOACCEL, ONE, ZERO + use global, only: cmfd, current_batch + use string, only: to_str + + 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 :: l ! iteration counter for leakage + real(8) :: flux ! temp variable for flux + real(8) :: leak1 ! group 1 leakage + real(8) :: leak2 ! group 2 leakage + + ! 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 + ZLOOP: do k = 1,nz + + YLOOP: do j = 1,ny + + XLOOP: do i = 1,nx + + ! Check for active mesh cell + if (allocated(cmfd%coremap)) then + if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) then + cycle + end if + end if + + ! Loop around energy groups + OUTGROUP: do h = 1,ng + + ! Get flux + flux = sum(cmfd % flux_rate(h,i,j,k,:)) + cmfd % flux(h,i,j,k) = flux + + ! Get total rr and convert to total xs + cmfd % totalxs(h,i,j,k) = sum(cmfd % total_rate(h,i,j,k,:)) / flux + + ! Get p1 scatter rr and convert to p1 scatter xs + cmfd % p1scattxs(h,i,j,k) = sum(cmfd % p1scatt_rate(h,i,j,k,:)) / flux + + ! Calculate diffusion coefficient + cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*(cmfd % totalxs(h,i,j,k) - & + cmfd % p1scattxs(h,i,j,k))) + + ! Begin loop to get energy out tallies + INGROUP: do g = 1, ng + + ! Get scattering + cmfd % scattxs(h,g,i,j,k) = sum(cmfd % scatt_rate(h,g,i,j,k,:)) / flux + + + ! Get nu-fission + cmfd % nfissxs(h,g,i,j,k) = sum(cmfd % nfiss_rate(h,g,i,j,k,:)) / flux + + end do INGROUP + + ! Surface currents + cmfd % current(1,h,i,j,k) = sum(cmfd % current_rate(1,h,i,j,k,:)) + cmfd % current(2,h,i,j,k) = sum(cmfd % current_rate(2,h,i,j,k,:)) + cmfd % current(3,h,i,j,k) = sum(cmfd % current_rate(3,h,i,j,k,:)) + cmfd % current(4,h,i,j,k) = sum(cmfd % current_rate(4,h,i,j,k,:)) + cmfd % current(5,h,i,j,k) = sum(cmfd % current_rate(5,h,i,j,k,:)) + cmfd % current(6,h,i,j,k) = sum(cmfd % current_rate(6,h,i,j,k,:)) + cmfd % current(7,h,i,j,k) = sum(cmfd % current_rate(7,h,i,j,k,:)) + cmfd % current(8,h,i,j,k) = sum(cmfd % current_rate(8,h,i,j,k,:)) + cmfd % current(9,h,i,j,k) = sum(cmfd % current_rate(9,h,i,j,k,:)) + cmfd % current(10,h,i,j,k) = sum(cmfd % current_rate(10,h,i,j,k,:)) + cmfd % current(11,h,i,j,k) = sum(cmfd % current_rate(11,h,i,j,k,:)) + cmfd % current(12,h,i,j,k) = sum(cmfd % current_rate(12,h,i,j,k,:)) + + end do OUTGROUP + + end do XLOOP + + end do YLOOP + + end do ZLOOP + #ifdef CMFD_DEBUG open(file='openmc_src_' // trim(to_str(current_batch)) // '.dat', unit=102) open(file='totalxs1_' // trim(to_str(current_batch)) // '.dat', unit=103) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index f89f8c3c46..87b4a89f86 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -97,8 +97,13 @@ contains ! If this is a restart run and we are just replaying batches leave if (restart_run .and. current_batch <= restart_batch) return - ! Check to reset tallies - if (cmfd_run .and. cmfd_reset % contains(current_batch)) then + ! Check to reset tallies with flusher + if (cmfd_run .and. (cmfd_reset % contains(current_batch))) then + call cmfd_tally_reset() + end if + + ! Check to reset tallies during CMFD with window + if (cmfd_run .and. cmfd_flush_every .and. current_batch > cmfd_begin) then call cmfd_tally_reset() end if diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 6a2dc8ed54..37352ac176 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -22,6 +22,15 @@ module cmfd_header ! Energy grid real(8), allocatable :: egrid(:) + ! Reaction rates + integer :: idx + real(8), allocatable :: flux_rate(:,:,:,:,:) + real(8), allocatable :: total_rate(:,:,:,:,:) + real(8), allocatable :: p1scatt_rate(:,:,:,:,:) + real(8), allocatable :: scatt_rate(:,:,:,:,:,:) + real(8), allocatable :: nfiss_rate(:,:,:,:,:,:) + real(8), allocatable :: current_rate(:,:,:,:,:,:) + ! Cross sections real(8), allocatable :: totalxs(:,:,:,:) real(8), allocatable :: p1scattxs(:,:,:,:) @@ -94,9 +103,10 @@ contains ! ALLOCATE_CMFD allocates all data in of cmfd type !============================================================================== - subroutine allocate_cmfd(this, n_batches) + subroutine allocate_cmfd(this, n_batches, n_save) integer, intent(in) :: n_batches ! number of batches in calc + integer, intent(in) :: n_save ! number of batches to save type(cmfd_type), intent(inout) :: this ! cmfd instance integer :: nx ! number of mesh cells in x direction @@ -104,12 +114,21 @@ contains integer :: nz ! number of mesh cells in z direction integer :: ng ! number of energy groups - ! Extract spatial and energy indices from object + ! Extract spatial and energy indices from object nx = this % indices(1) ny = this % indices(2) nz = this % indices(3) ng = this % indices(4) + ! Allocate all rates + if (.not. allocated(this % flux_rate)) allocate(this % flux_rate(ng,nx,ny,nz,n_save)) + if (.not. allocated(this % total_rate)) allocate(this % total_rate(ng,nx,ny,nz,n_save)) + if (.not. allocated(this % p1scatt_rate)) allocate(this % p1scatt_rate(ng,nx,ny,nz,n_save)) + if (.not. allocated(this % scatt_rate)) allocate(this % scatt_rate(ng,ng,nx,ny,nz,n_save)) + if (.not. allocated(this % nfiss_rate)) allocate(this % nfiss_rate(ng,ng,nx,ny,nz,n_save)) + if (.not. allocated(this % current_rate)) allocate(this % current_rate(12,ng,nx,ny,nz,n_save)) + + ! Allocate flux, cross sections and diffusion coefficient if (.not. allocated(this % flux)) allocate(this % flux(ng,nx,ny,nz)) if (.not. allocated(this % totalxs)) allocate(this % totalxs(ng,nx,ny,nz)) @@ -149,11 +168,17 @@ contains this % p1scattxs = ZERO this % scattxs = ZERO this % nfissxs = ZERO + this % flux_rate = ZERO + this % total_rate = ZERO + this % p1scatt_rate = ZERO + this % scatt_rate = ZERO + this % nfiss_rate = ZERO this % diffcof = ZERO this % dtilde = ZERO this % dhat = ZERO this % hxyz = ZERO this % current = ZERO + this % current_rate = ZERO this % cmfd_src = ZERO this % openmc_src = ZERO this % sourcecounts = ZERO @@ -164,6 +189,9 @@ contains this % k_cmfd = ZERO this % entropy = ZERO + ! Set starting index to 1 + this % idx = 1 + end subroutine allocate_cmfd !=============================================================================== @@ -182,6 +210,12 @@ contains if (allocated(this % diffcof)) deallocate(this % diffcof) if (allocated(this % current)) deallocate(this % current) if (allocated(this % flux)) deallocate(this % flux) + if (allocated(this % total_rate)) deallocate(this % total_rate) + if (allocated(this % p1scatt_rate)) deallocate(this % p1scatt_rate) + if (allocated(this % scatt_rate)) deallocate(this % scatt_rate) + if (allocated(this % nfiss_rate)) deallocate(this % nfiss_rate) + if (allocated(this % current_rate)) deallocate(this % current_rate) + if (allocated(this % flux_rate)) deallocate(this % flux_rate) if (allocated(this % dtilde)) deallocate(this % dtilde) if (allocated(this % dhat)) deallocate(this % dhat) if (allocated(this % hxyz)) deallocate(this % hxyz) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 4ef30f231f..9bfa5df4de 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -54,7 +54,7 @@ contains call time_cmfdsolve % reset() ! Allocate cmfd object - call allocate_cmfd(cmfd, n_batches) + call allocate_cmfd(cmfd, n_batches, cmfd_n_save) end subroutine configure_cmfd @@ -295,6 +295,12 @@ contains end select end if + ! Check for saving batches + if (check_for_node(doc, "n_save")) then + call get_node_value(doc, "n_save", cmfd_n_save) + cmfd_flush_every = .true. + end if + ! Create tally objects call create_cmfd_tally(doc, estimator) diff --git a/src/global.F90 b/src/global.F90 index 1e8526c475..9e9096e380 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -371,6 +371,10 @@ module global real(8) :: cmfd_atoli = 1.e-10_8 real(8) :: cmfd_rtoli = 1.e-5_8 + ! Number of batches to save for in CMFD + integer :: cmfd_n_save = 1 + logical :: cmfd_flush_every = .false. + ! Information about state points to be written integer :: n_state_points = 0 type(SetInt) :: statepoint_batch