From a6fa726e18daef8fe0daf70ee7ff73e7645da4ce Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 28 Feb 2014 23:26:02 -0500 Subject: [PATCH 01/25] added option to write out initial source guess --- src/global.F90 | 3 +++ src/input_xml.F90 | 8 ++++++++ src/source.F90 | 11 +++++++++++ 3 files changed, 22 insertions(+) diff --git a/src/global.F90 b/src/global.F90 index b07132270d..7ec2f0f19c 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -299,6 +299,9 @@ module global ! Particle restart run logical :: particle_restart_run = .false. + ! Write out initial source + logical :: write_initial_source = .false. + ! ============================================================================ ! CMFD VARIABLES diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 039e0419c3..842c104b36 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -262,6 +262,14 @@ contains call fatal_error() end if + ! Check if we want to write out source + if (check_for_node(node_source, "write_initial")) then + call get_node_value(node_source, "write_initial", temp_str) + call lower_case(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + write_initial_source = .true. + end if + ! Check for external source file if (check_for_node(node_source, "file")) then ! Copy path of source file diff --git a/src/source.F90 b/src/source.F90 index 9dcfd01445..090ea93d1e 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -27,6 +27,7 @@ contains subroutine initialize_source() + character(MAX_FILE_LEN) :: filename integer(8) :: i ! loop index over bank sites integer(8) :: id ! particle id integer(4) :: itmp ! temporary integer @@ -76,6 +77,16 @@ contains end do end if + ! Write out initial source + if (write_initial_source) then + message = 'Writing out initial source guess...' + call write_message(1) + filename = 'initial_source.h5' + call sp % file_create(filename, serial = .false.) + call sp % write_source_bank() + call sp % file_close() + end if + end subroutine initialize_source !=============================================================================== From 0a31e55e95539ad3b6fee25e16c6586ca55891f4 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Sun, 2 Mar 2014 18:28:14 -0500 Subject: [PATCH 02/25] flushing now performed with sets --- src/cmfd_execute.F90 | 20 ++++---------------- src/cmfd_input.F90 | 24 +++++++++++++++--------- src/global.F90 | 8 +++----- 3 files changed, 22 insertions(+), 30 deletions(-) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 35485a8945..e4e849436f 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -82,8 +82,8 @@ contains subroutine cmfd_init_batch() - use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & - cmfd_inact_flush, cmfd_act_flush, cmfd_run, & + use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & + cmfd_reset, cmfd_run, & current_batch, cmfd_hold_weights ! Check to activate CMFD diffusion and possible feedback @@ -97,21 +97,9 @@ 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 flush cmfd tallies for active batches, no more inactive flush - if (cmfd_run .and. cmfd_act_flush == current_batch) then + ! Check to reset tallies + if (cmfd_run .and. cmfd_reset % contains(current_batch)) then call cmfd_tally_reset() - cmfd_tally_on = .true. - cmfd_inact_flush(2) = -1 - end if - - ! Check to flush cmfd tallies during inactive batches (>= on number of - ! flushes important as the code will flush on the first batch which we - ! dont want to count) - if (cmfd_run .and. mod(current_batch,cmfd_inact_flush(1)) & - == 0 .and. cmfd_inact_flush(2) > 0 .and. cmfd_begin < current_batch) then - cmfd_hold_weights = .true. - call cmfd_tally_reset() - cmfd_inact_flush(2) = cmfd_inact_flush(2) - 1 end if end subroutine cmfd_init_batch diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 61d6f6eaba..93454ce4c2 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -68,6 +68,7 @@ contains integer :: ng integer, allocatable :: iarray(:) + integer, allocatable :: int_array(:) logical :: file_exists ! does cmfd.xml exist? logical :: found character(MAX_LINE_LEN) :: filename @@ -216,15 +217,20 @@ contains cmfd_tally_on = .false. end if - ! Inactive batch flush window - if (check_for_node(doc, "inactive_flush")) & - call get_node_value(doc, "inactive_flush", cmfd_inact_flush(1)) - if (check_for_node(doc, "num_flushes")) & - call get_node_value(doc, "num_flushes", cmfd_inact_flush(2)) - - ! Last flush before active batches - if (check_for_node(doc, "active_flush")) & - call get_node_value(doc, "active_flush", cmfd_act_flush) + ! Check for cmfd tally resets + if (check_for_node(doc, "tally_reset")) then + n_cmfd_resets = get_arraysize_integer(doc, "tally_reset") + else + n_cmfd_resets = 0 + end if + if (n_cmfd_resets > 0) then + allocate(int_array(n_cmfd_resets)) + call get_node_array(doc, "tally_reset", int_array) + do i = 1, n_cmfd_resets + call cmfd_reset % add(int_array(i)) + end do + deallocate(int_array) + end if ! Get display if (check_for_node(doc, "display")) & diff --git a/src/global.F90 b/src/global.F90 index 7ec2f0f19c..02a87051f4 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -347,11 +347,9 @@ module global ! Batch to begin cmfd integer :: cmfd_begin = 1 - ! When and how long to flush cmfd tallies during inactive batches - integer :: cmfd_inact_flush(2) = (/9999,1/) - - ! Batch to last flush before active batches - integer :: cmfd_act_flush = 0 + ! Tally reset list + integer :: n_cmfd_resets + type(SetInt) :: cmfd_reset ! Compute effective downscatter cross section logical :: cmfd_downscatter = .false. From 16a20d10e816b4b75f03e2d1df40d9a4f1955ece Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Sun, 2 Mar 2014 18:52:23 -0500 Subject: [PATCH 03/25] added flag for cmfd debug printing --- src/CMakeLists.txt | 10 +++++ src/cmfd_data.F90 | 90 +++++++++++++++++++++++++++++++++++++++++++- src/cmfd_execute.F90 | 25 +++++++++++- 3 files changed, 123 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fc6c6630cc..033ce2a020 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -26,6 +26,7 @@ option(debug "Compile with debug flags" OFF) option(optimize "Turn on all compiler optimization flags" OFF) option(verbose "Create verbose Makefiles" OFF) option(coverage "Compile with flags" OFF) +option(cmfd_debug "Compile with cmfd debug flags" OFF) if (verbose) set(CMAKE_VERBOSE_MAKEFILE on) @@ -214,6 +215,15 @@ if(GIT_SHA1_SUCCESS EQUAL 0) add_definitions(-DGIT_SHA1="${GIT_SHA1}") endif() +#=============================================================================== +# CMFD debug flags +#=============================================================================== + +if (cmfd_debug) + message("-- CMFD Debug flags on") + add_definitions(-DCMFD_DEBUG) +endif() + #=============================================================================== # FoX Fortran XML Library #=============================================================================== diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index b5612dcd05..7962760cbf 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -56,7 +56,7 @@ contains ONE, TINY_BIT use error, only: fatal_error use global, only: cmfd, message, n_cmfd_tallies, cmfd_tallies, meshes,& - matching_bins + matching_bins, current_batch use mesh, only: mesh_indices_to_bin use mesh_header, only: StructuredMesh use string, only: to_str @@ -71,6 +71,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 :: ital ! tally object index integer :: ijk(3) ! indices for mesh cell integer :: score_index ! index to pull from tally object @@ -80,6 +81,8 @@ contains 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 @@ -301,6 +304,91 @@ contains if (associated(t)) nullify(t) if (associated(m)) nullify(m) +#ifdef CMFD_DEBUG + open(file='cmfd_src_' // trim(to_str(current_batch)) // '.dat', unit=100) + open(file='source_bank_' // trim(to_str(current_batch)) // '.dat', unit=101) + open(file='openmc_src_' // trim(to_str(current_batch)) // '.dat', unit=102) + open(file='totalxs1_' // trim(to_str(current_batch)) // '.dat', unit=103) + open(file='totalxs2_' // trim(to_str(current_batch)) // '.dat', unit=104) + open(file='p1scattxs1_' // trim(to_str(current_batch)) // '.dat', unit=105) + open(file='p1scattxs2_' // trim(to_str(current_batch)) // '.dat', unit=106) + open(file='scattxs11_' // trim(to_str(current_batch)) // '.dat', unit=107) + open(file='scattxs12_' // trim(to_str(current_batch)) // '.dat', unit=108) + open(file='scattxs21_' // trim(to_str(current_batch)) // '.dat', unit=109) + open(file='scattxs22_' // trim(to_str(current_batch)) // '.dat', unit=110) + open(file='nufissxs11_' // trim(to_str(current_batch)) // '.dat', unit=111) + open(file='nufissxs12_' // trim(to_str(current_batch)) // '.dat', unit=112) + open(file='nufissxs21_' // trim(to_str(current_batch)) // '.dat', unit=113) + open(file='nufissxs22_' // trim(to_str(current_batch)) // '.dat', unit=114) + open(file='diff_coef1_' // trim(to_str(current_batch)) // '.dat', unit=115) + open(file='diff_coef2_' // trim(to_str(current_batch)) // '.dat', unit=116) + open(file='flux1_' // trim(to_str(current_batch)) // '.dat', unit=117) + open(file='flux2_' // trim(to_str(current_batch)) // '.dat', unit=118) + open(file='leak1_' // trim(to_str(current_batch)) // '.dat', unit=119) + open(file='leak2_' // trim(to_str(current_batch)) // '.dat', unit=120) + + do i = 1, nx + do j = 1,ny + leak1 = ZERO + leak2 = ZERO + do l = 1, 3 + leak1 = leak1 + ((cmfd % current(4*l,1,i,j,1) - & + cmfd % current(4*l-1,1,i,j,1))) - & + ((cmfd % current(4*l-2,1,i,j,1) - & + cmfd % current(4*l-3,1,i,j,1))) + end do + do l = 1, 3 + leak2 = leak2 + ((cmfd % current(4*l,2,i,j,1) - & + cmfd % current(4*l-1,2,i,j,1))) - & + ((cmfd % current(4*l-2,2,i,j,1) - & + cmfd % current(4*l-3,2,i,j,1))) + end do + write(100,*) cmfd % cmfd_src(1,i,j,1) + write(101,*) cmfd % sourcecounts(1,i,j,1) + write(102,*) cmfd % openmc_src(1,i,j,1) + write(103,*) cmfd % totalxs(1,i,j,1) + write(104,*) cmfd % totalxs(2,i,j,1) + write(105,*) cmfd % p1scattxs(1,i,j,1) + write(106,*) cmfd % p1scattxs(2,i,j,1) + write(107,*) cmfd % scattxs(1,1,i,j,1) + write(108,*) cmfd % scattxs(1,2,i,j,1) + write(109,*) cmfd % scattxs(2,1,i,j,1) + write(110,*) cmfd % scattxs(2,2,i,j,1) + write(111,*) cmfd % nfissxs(1,1,i,j,1) + write(112,*) cmfd % nfissxs(1,2,i,j,1) + write(113,*) cmfd % nfissxs(2,1,i,j,1) + write(114,*) cmfd % nfissxs(2,2,i,j,1) + write(115,*) cmfd % diffcof(1,i,j,1) + write(116,*) cmfd % diffcof(2,i,j,1) + write(117,*) cmfd % flux(1,i,j,1) + write(118,*) cmfd % flux(2,i,j,1) + write(119,*) leak1 + write(120,*) leak2 + end do + end do + close(100) + close(101) + close(102) + close(103) + close(104) + close(105) + close(106) + close(107) + close(108) + close(109) + close(110) + close(111) + close(112) + close(113) + close(114) + close(115) + close(116) + close(117) + close(118) + close(119) + close(120) +#endif + end subroutine compute_xs !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index e4e849436f..a1b8853042 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -69,7 +69,7 @@ contains call calc_fission_source() ! calculate weight factors - if (cmfd_feedback) call cmfd_reweight(.true.) + call cmfd_reweight(.true.) ! stop cmfd timer if (master) call time_cmfd % stop() @@ -132,6 +132,7 @@ contains use constants, only: CMFD_NOACCEL, ZERO, TWO use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -238,6 +239,15 @@ contains cmfd % src_cmp(current_batch) = sqrt(ONE/cmfd % norm * & sum((cmfd % cmfd_src - cmfd % openmc_src)**2)) +#ifdef CMFD_DEBUG + open(file='cmfd_src_' // trim(to_str(current_batch)) // '.dat', unit=100) + do i = 1, nx + do j = 1,ny + write(100,*) cmfd % sourcecounts(1,i,j,1) + end do + end do + close(100) +#endif end if #ifdef MPI @@ -260,6 +270,7 @@ contains use mesh_header, only: StructuredMesh use mesh, only: count_bank_sites, get_mesh_indices use search, only: binary_search + use string, only: to_str #ifdef MPI use global, only: mpi_err @@ -273,6 +284,7 @@ contains integer :: nz ! maximum number of cells in z direction integer :: ng ! maximum number of energy groups integer :: i ! iteration counter + integer :: j ! iteration counter integer :: ijk(3) ! spatial bin location integer :: e_bin ! energy bin of source particle integer :: n_groups ! number of energy groups @@ -327,8 +339,19 @@ contains cmfd % weightfactors = cmfd % cmfd_src/sum(cmfd % cmfd_src)* & sum(cmfd % sourcecounts) / cmfd % sourcecounts end where +#ifdef CMFD_DEBUG + open(file='source_bank_' // trim(to_str(current_batch)) // '.dat', unit=101) + do i = 1, nx + do j = 1,ny + write(101,*) cmfd % sourcecounts(1,i,j,1) + end do + end do + close(101) +#endif end if +if (.not. cmfd_feedback) return + ! Broadcast weight factors to all procs #ifdef MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & From cf1ace33032acb370b383c6e72bf7b33b830ef3c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Sun, 2 Mar 2014 19:39:16 -0500 Subject: [PATCH 04/25] added option to reset dhats --- src/cmfd_data.F90 | 6 ------ src/cmfd_execute.F90 | 2 +- src/cmfd_input.F90 | 8 ++++++++ 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 7962760cbf..df6b38e910 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -305,8 +305,6 @@ contains if (associated(m)) nullify(m) #ifdef CMFD_DEBUG - open(file='cmfd_src_' // trim(to_str(current_batch)) // '.dat', unit=100) - open(file='source_bank_' // trim(to_str(current_batch)) // '.dat', unit=101) open(file='openmc_src_' // trim(to_str(current_batch)) // '.dat', unit=102) open(file='totalxs1_' // trim(to_str(current_batch)) // '.dat', unit=103) open(file='totalxs2_' // trim(to_str(current_batch)) // '.dat', unit=104) @@ -343,8 +341,6 @@ contains ((cmfd % current(4*l-2,2,i,j,1) - & cmfd % current(4*l-3,2,i,j,1))) end do - write(100,*) cmfd % cmfd_src(1,i,j,1) - write(101,*) cmfd % sourcecounts(1,i,j,1) write(102,*) cmfd % openmc_src(1,i,j,1) write(103,*) cmfd % totalxs(1,i,j,1) write(104,*) cmfd % totalxs(2,i,j,1) @@ -366,8 +362,6 @@ contains write(120,*) leak2 end do end do - close(100) - close(101) close(102) close(103) close(104) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index a1b8853042..f89f8c3c46 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -243,7 +243,7 @@ contains open(file='cmfd_src_' // trim(to_str(current_batch)) // '.dat', unit=100) do i = 1, nx do j = 1,ny - write(100,*) cmfd % sourcecounts(1,i,j,1) + write(100,*) cmfd % cmfd_src(1,i,j,1) end do end do close(100) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 93454ce4c2..39c780173f 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -165,6 +165,14 @@ contains cmfd_downscatter = .true. end if + ! Run an adjoint calc + if (check_for_node(doc, "dhat_reset")) then + call get_node_value(doc, "dhat_reset", temp_str) + call lower_case(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + dhat_reset = .true. + end if + ! Set the solver type if (check_for_node(doc, "solver")) & call get_node_value(doc, "solver", cmfd_solver_type) From ccbd5b712817b09f4f8ec41beaaa3e18ba87cb23 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Sun, 2 Mar 2014 19:42:15 -0500 Subject: [PATCH 05/25] added fix neutron balance option --- src/cmfd_input.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 39c780173f..2cd3d1f619 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -165,7 +165,7 @@ contains cmfd_downscatter = .true. end if - ! Run an adjoint calc + ! Reset dhat parameters if (check_for_node(doc, "dhat_reset")) then call get_node_value(doc, "dhat_reset", temp_str) call lower_case(temp_str) @@ -173,6 +173,14 @@ contains dhat_reset = .true. end if + ! Enforce perfect neutorn balance + if (check_for_node(doc, "balance")) then + call get_node_value(doc, "balance", temp_str) + call lower_case(temp_str) + if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & + cmfd_fix_balance = .true. + end if + ! Set the solver type if (check_for_node(doc, "solver")) & call get_node_value(doc, "solver", cmfd_solver_type) From 4de45e7b08b9d6a71a8d1c38b9de2bdd066f4809 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 3 Mar 2014 11:16:50 -0500 Subject: [PATCH 06/25] dhats can be reset from cmfd input --- src/cmfd_data.F90 | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index df6b38e910..1a35ac09bc 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -10,8 +10,6 @@ module cmfd_data private public :: set_up_cmfd, neutron_balance - logical :: dhat_reset = .false. - contains !============================================================================== @@ -705,7 +703,8 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use global, only: cmfd, cmfd_coremap + use global, only: cmfd, cmfd_coremap, message, dhat_reset + use output, only: write_message integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -825,7 +824,9 @@ contains cmfd%dhat(l,g,i,j,k) = dhat ! check for dhat reset - if (dhat_reset) cmfd%dhat(l,g,i,j,k) = ZERO + if (dhat_reset) then + cmfd%dhat(l,g,i,j,k) = ZERO + end if end do LEAK @@ -837,6 +838,12 @@ contains end do ZLOOP + ! write that dhats are zero + if (dhat_reset) then + message = 'Dhats reset to zero.' + call write_message(1) + end if + end subroutine compute_dhat !=============================================================================== @@ -914,6 +921,16 @@ contains real(8) :: siga1 ! group 1 abs xs real(8) :: siga2 ! group 2 abs xs real(8) :: sigs12_eff ! effective downscatter xs + real(8) :: a ! matrix (1,1) element for balance + real(8) :: b ! matrix (1,2) element for balance + real(8) :: c ! matrix (2,1) element for balance + real(8) :: d ! matrix (2,2) element for balance + real(8) :: r1 ! right hand side element 1 + real(8) :: r2 ! right hand side element 2 + real(8) :: det ! determinant of balance matrix + + message = 'Correcting neutron balance' + call write_message(1) ! Extract spatial and energy indices from object nx = cmfd % indices(1) From d1e46a22ec3cd49c068dc96ab1ab367a47f9323e Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Sat, 15 Mar 2014 08:18:07 -0400 Subject: [PATCH 07/25] added keff balance --- src/cmfd_data.F90 | 52 +++++++++++++++++++++++++++++++++------------ src/cmfd_header.F90 | 3 +++ 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 1a35ac09bc..dc5db8352c 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -31,6 +31,9 @@ contains ! Compute effective downscatter cross section if (cmfd_downscatter) call compute_effective_downscatter() + ! Compute perfect balance + if (cmfd_fix_balance) call fix_neutron_balance() + ! Check neutron balance call neutron_balance() @@ -104,6 +107,8 @@ contains cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width + cmfd % keff_bal = ZERO + ! Begin loop around tallies TAL: do ital = 1, n_cmfd_tallies @@ -212,6 +217,7 @@ contains ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & t % results(2,score_index) % sum + cmfd % keff_bal = cmfd % keff_bal + t % results(2,score_index) % sum / dble(t % n_realizations) end do INGROUP @@ -703,8 +709,9 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use global, only: cmfd, cmfd_coremap, message, dhat_reset + use global, only: cmfd, cmfd_coremap, message, dhat_reset, current_batch use output, only: write_message + use string, only: to_str integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -738,6 +745,9 @@ contains nxyz(1,:) = (/1,nx/) nxyz(2,:) = (/1,ny/) nxyz(3,:) = (/1,nz/) +#ifdef CMFD_DEBUG + open(file='cmfd_dhat_' // trim(to_str(current_batch)) // '.dat', unit=125) +#endif ! Geting loop over group and spatial indices ZLOOP: do k = 1,nz @@ -828,6 +838,9 @@ contains cmfd%dhat(l,g,i,j,k) = ZERO end if +#ifdef CMFD_DEBUG + write(125,*) dhat +#endif end do LEAK end do GROUP @@ -844,6 +857,10 @@ contains call write_message(1) end if +#ifdef CMFD_DEBUG + close(unit=125) +#endif + end subroutine compute_dhat !=============================================================================== @@ -889,12 +906,12 @@ contains !=============================================================================== ! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance !=============================================================================== -#ifdef DEVELOPMENTAL + subroutine fix_neutron_balance() use constants, only: ONE, ZERO, CMFD_NOACCEL - use global, only: cmfd, keff - use, intrinsic :: ISO_FORTRAN_ENV + use global, only: cmfd, message + use output, only: write_message integer :: nx ! number of mesh cells in x direction integer :: ny ! number of mesh cells in y direction @@ -928,6 +945,7 @@ contains real(8) :: r1 ! right hand side element 1 real(8) :: r2 ! right hand side element 2 real(8) :: det ! determinant of balance matrix + real(8) :: keff_bal ! keffective for balance eq. message = 'Correcting neutron balance' call write_message(1) @@ -938,6 +956,10 @@ contains nz = cmfd % indices(3) ng = cmfd % indices(4) + keff_bal = cmfd % keff_bal + + print *, 'KEFF BALANCE is:', keff_bal + ! Return if not two groups if (ng /= 2) return @@ -995,14 +1017,16 @@ contains siga1 = sigt1 - sigs11 - sigs12 siga2 = sigt2 - sigs22 - sigs21 - ! Compute effective downscatter xs - sigs12_eff = (ONE/keff*nsigf11*flux1 - leak1 - siga1*flux1 & - - ONE/keff*nsigf21/siga2*leak2 ) / ( flux1*(ONE & - - ONE/keff*nsigf21/siga2)) - - ! Redefine flux 2 - flux2 = (sigs12_eff*flux1 - leak2)/siga2 - cmfd % flux(2,i,j,k) = flux2 + ! Create matrix and solve for effective downscatter and thermal flux + a = flux1 + b = -ONE/keff_bal*nsigf21 + c = flux1 + d = -siga2 + ONE/keff_bal*nsigf22 + det = a*d - b*c + r1 = ONE/keff_bal*nsigf11*flux1 - siga1*flux1 - leak1 + r2 = leak2 - ONE/keff_bal*nsigf12*flux1 + sigs12_eff = ONE/det*(d*r1 - b*r2) + flux2 = ONE/det*(a*r2 - c*r1) ! Recompute total cross sections (use effective and no upscattering) sigt1 = siga1 + sigs11 + sigs12_eff @@ -1018,6 +1042,9 @@ contains ! Zero out upscatter cross section cmfd % scattxs(2,1,i,j,k) = ZERO + ! Record thermal flux + cmfd % flux(2,i,j,k) = flux2 + end do XLOOP end do YLOOP @@ -1025,7 +1052,6 @@ contains end do ZLOOP end subroutine fix_neutron_balance -#endif !=============================================================================== ! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 73bc0e0163..c9bc0f32af 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -83,6 +83,9 @@ module cmfd_header ! List of CMFD k real(8), allocatable :: k_cmfd(:) + ! Balance keff + real(8) :: keff_bal + end type cmfd_type contains From d26b2da8bc4936a0e3d2d50c74205f9e4f6e6fad Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 11:51:52 -0400 Subject: [PATCH 08/25] updated documentation and removed some unncessary vars --- docs/source/usersguide/input.rst | 136 +++++++++++++++++++++---------- src/cmfd_execute.F90 | 3 +- src/cmfd_input.F90 | 14 ++-- src/global.F90 | 1 - 4 files changed, 97 insertions(+), 57 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c9aa8cc866..301f7f52de 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1238,14 +1238,26 @@ Currently, it allows users to accelerate fission source convergence during inactive neutron batches. To run CMFD, the ```` element in ``settings.xml`` should be set to "true". -```` Element --------------------------- +```` Element +-------------------- -The ```` element controls the batch where CMFD tallies should be -reset. CMFD tallies should be reset before active batches so they are accumulated -without bias. +The ```` element specifies the absolute inner tolerance on the Gauss-Seidel +when performing CMFD calculations. It is only used in the standalone CMFD solver +and not when PETSc is active. + + *Default*: 1.e-10 + +```` Element +------------------- + +The ```` element controls whether exact neutron balance should be enforced +from CMFD tallies before creating CMFD matrices. This changes effective downscatter +cross section and thermal flux to create exact balance. It should be noted that +this option has led to instabilities when performing CMFD. It can be turned on +with "true" and off with "false". + + *Default*: false - *Default*: 0 ```` Element ------------------- @@ -1268,7 +1280,25 @@ The ```` element sets one additional CMFD output column. Options are: * "source" - prints the RMS [%] between the OpenMC fission source and CMFD fission source. - *Default*: None + *Default*: balance + +```` Element +------------------- + +The ```` element controls whether dhat nonlinear CMFD parameters +should be reset to zero before solving CMFD eigenproblem. It can be turned on with +"true" and off with "false". + + *Default*: false + +```` Element +------------------- + +The ```` element controls whether an effective downscatter cross +section should be used when using 2-group CMFD. It can be turned on with "true" +and off with "false". + + *Default*: false ```` Element ---------------------- @@ -1279,25 +1309,6 @@ It can be turned on with "true" and off with "false". *Default*: false -```` Element ----------------------- - -The ```` element controls if cmfd tallies should be accumulated -during inactive batches. For some applications, CMFD tallies may not be -needed until the start of active batches. This option can be turned on -with "true" and off with "false" - - *Default*: true - -```` Element ----------------------------- - -The ```` element controls when CMFD tallies are reset during -inactive batches. The integer set here is the interval at which this reset -occurs. The amout of resets is controlled with the ```` element. - - *Defualt*: 9999 - ```` Element ------------------------- @@ -1305,9 +1316,16 @@ The ```` element is used to view the convergence of linear GMRES iterations in PETSc. This option can be turned on with "true" and turned off with "false". - *Default*: false +```` Element +-------------------- + +The ```` element specifies the tolerance on the eigenvalue when performing +CMFD power iteration. + + *Default*: 1.e-8 + ```` Element ------------------ @@ -1376,14 +1394,6 @@ not impact the calculation. *Default*: 1.0 -```` Element -------------------------- - -The ```` element controls the number of CMFD tally resets that -occur during inactive CMFD batches. - - *Default*: 9999 - ```` Element --------------------------- @@ -1392,20 +1402,21 @@ This option can be turned on with "true" and turned off with "false". *Default*: false +```` Element +-------------------- + +The ```` element specifies the relative inner tolerance on the Gauss-Seidel +when performing CMFD calculations. It is only used in the standalone CMFD solver +and not when PETSc is active. + + *Default*: 1.e-5 + ```` Element ------------------------- The ```` element can be turned on with "true" to have an adjoint -calculation be performed on the last batch when CMFD is active. - - *Default*: false - -```` Element --------------------------- - -The ```` element is used to view the convergence of the nonlinear SNES -function in PETSc. This option can be turned on with "true" and turned off with "false". - +calculation be performed on the last batch when CMFD is active. OpenMC should be +compiled with PETSc when using this option. *Default*: false @@ -1418,6 +1429,41 @@ By setting "power", power iteration is used and by setting "jfnk", JFNK is used. *Default*: power +```` Element +-------------------- + +The ```` element specifies an optional Wielandt shift parameter for +accelerating power iterations. It can only be used when PETSc is not active. +It is by default very large so the impact of the shift is effectively zero. + + *Default*: 1e6 + +```` Element +-------------------- + +The ```` element specifies an optional spectral radius that can be set to +accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration +solve. Note this is only used in the standalone CMFD solver and does not affect +the calculation when PETSc is active. + + *Default*: power + +```` Element +-------------------- + +The ```` element specifies the tolerance on the fission source when performing +CMFD power iteration. + + *Default*: 1.e-8 + +```` Element +-------------------- + +The ```` element contains a list of batch numbers in which CMFD tallies +should be reset. + + *Default*: None + ```` Element ---------------------------- diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index f89f8c3c46..fe19bd95b7 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -82,7 +82,7 @@ contains subroutine cmfd_init_batch() - use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & + use global, only: cmfd_begin, cmfd_on, & cmfd_reset, cmfd_run, & current_batch, cmfd_hold_weights @@ -91,7 +91,6 @@ contains ! accumulated if (cmfd_run .and. cmfd_begin == current_batch) then cmfd_on = .true. - cmfd_tally_on = .true. end if ! If this is a restart run and we are just replaying batches leave diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 2cd3d1f619..0523325c9a 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -173,7 +173,7 @@ contains dhat_reset = .true. end if - ! Enforce perfect neutorn balance + ! Enforce perfect neutron balance if (check_for_node(doc, "balance")) then call get_node_value(doc, "balance", temp_str) call lower_case(temp_str) @@ -218,6 +218,10 @@ contains call get_node_value(doc, "run_adjoint", temp_str) call lower_case(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & +#ifndef PETSC + message = 'Must use PETSc when running adjoint option.' + call fatal_error() +#endif cmfd_run_adjoint = .true. end if @@ -225,14 +229,6 @@ contains if (check_for_node(doc, "begin")) & call get_node_value(doc, "begin", cmfd_begin) - ! Tally during inactive batches - if (check_for_node(doc, "inactive")) then - call get_node_value(doc, "inactive", temp_str) - call lower_case(temp_str) - if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & - cmfd_tally_on = .false. - end if - ! Check for cmfd tally resets if (check_for_node(doc, "tally_reset")) then n_cmfd_resets = get_arraysize_integer(doc, "tally_reset") diff --git a/src/global.F90 b/src/global.F90 index 02a87051f4..3b9b8bfeea 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -367,7 +367,6 @@ module global ! CMFD run logicals logical :: cmfd_on = .false. - logical :: cmfd_tally_on = .true. ! CMFD display info character(len=25) :: cmfd_display = 'balance' From 445f3910f8cf06d5f81598164bdea5668270a650 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 9 Sep 2014 11:56:21 -0400 Subject: [PATCH 09/25] fixed titles in cmfd section --- docs/source/usersguide/input.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 301f7f52de..cc2d5e7597 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1248,7 +1248,7 @@ and not when PETSc is active. *Default*: 1.e-10 ```` Element -------------------- +--------------------- The ```` element controls whether exact neutron balance should be enforced from CMFD tallies before creating CMFD matrices. This changes effective downscatter @@ -1283,7 +1283,7 @@ The ```` element sets one additional CMFD output column. Options are: *Default*: balance ```` Element -------------------- +------------------------ The ```` element controls whether dhat nonlinear CMFD parameters should be reset to zero before solving CMFD eigenproblem. It can be turned on with @@ -1292,7 +1292,7 @@ should be reset to zero before solving CMFD eigenproblem. It can be turned on wi *Default*: false ```` Element -------------------- +------------------------- The ```` element controls whether an effective downscatter cross section should be used when using 2-group CMFD. It can be turned on with "true" @@ -1439,7 +1439,7 @@ It is by default very large so the impact of the shift is effectively zero. *Default*: 1e6 ```` Element --------------------- +---------------------- The ```` element specifies an optional spectral radius that can be set to accelerate the convergence of Gauss-Seidel iterations during CMFD power iteration @@ -1449,7 +1449,7 @@ the calculation when PETSc is active. *Default*: power ```` Element --------------------- +------------------ The ```` element specifies the tolerance on the fission source when performing CMFD power iteration. @@ -1457,7 +1457,7 @@ CMFD power iteration. *Default*: 1.e-8 ```` Element --------------------- +------------------------- The ```` element contains a list of batch numbers in which CMFD tallies should be reset. From 2c98fd3d6e486c916dc099a1ebfaafced0199940 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 15 Apr 2014 17:26:45 -0400 Subject: [PATCH 10/25] added capability for more than 2 groups --- src/cmfd_solver.F90 | 105 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 99 insertions(+), 6 deletions(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 11989ff4f6..dd343d87fc 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -132,7 +132,7 @@ contains guess = ONE phi_n % val = guess phi_o % val = guess - k_n = keff + k_n = keff k_o = k_n dw = cmfd_shift k_s = k_o + dw @@ -162,10 +162,9 @@ contains case(1) cmfd_linsolver => cmfd_linsolver_1g case(2) - cmfd_linsolver => cmfd_linsolver_2g + cmfd_linsolver => cmfd_linsolver_ng case default - message = 'Must use PETSc for more than 2 groups' - call fatal_error() + cmfd_linsolver => cmfd_linsolver_ng end select ! Set tolerances @@ -306,7 +305,7 @@ contains subroutine convergence(iter, innerits) - use constants, only: ONE, TINY_BIT + use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV @@ -320,7 +319,7 @@ contains kerr = abs(k_o - k_n)/k_n ! Calculate max error in source - where (s_n % val > TINY_BIT) + where (s_n % val > ZERO) serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 end where serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) @@ -582,6 +581,100 @@ contains end subroutine cmfd_linsolver_2g +!=============================================================================== +! CMFD_LINSOLVER_ng solves the CMFD linear system +!=============================================================================== + + subroutine cmfd_linsolver_ng(A, b, x, tol, its) + + use constants, only: ONE, ZERO + use global, only: cmfd, cmfd_spectral + + type(Matrix) :: A ! coefficient matrix + type(Vector) :: b ! right hand side vector + type(Vector) :: x ! unknown vector + real(8) :: tol ! tolerance on final error + integer :: its ! number of inner iterations + + integer :: g ! group index + integer :: i ! loop counter for x + integer :: j ! loop counter for y + integer :: k ! loop counter for z + integer :: n ! total size of vector + integer :: nx ! maximum dimension in x direction + integer :: ny ! maximum dimension in y direction + integer :: nz ! maximum dimension in z direction + integer :: ng ! number of energy groups + integer :: igs ! Gauss-Seidel iteration counter + integer :: irow ! row iteration + integer :: icol ! iteration counter over columns + integer :: didx ! index for diagonal component + logical :: found ! did we find col + real(8) :: tmp1 ! temporary sum g1 + real(8) :: x1 ! new g1 value of x + real(8) :: err ! error in convergence of solution + real(8) :: w ! overrelaxation parameter + type(Vector) :: tmpx ! temporary solution vector + + ! Set overrelaxation parameter + w = ONE + + ! Dimensions + ng = 1 + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + n = A % n + + ! Perform Gauss Seidel iterations + GS: do igs = 1, 10000 + + ! Copy over x vector + call tmpx % copy(x) + + ! Begin loop around matrix rows + ROWS: do irow = 1, n + + ! Get spatial location + call matrix_to_indices(irow, g, i, j, k, ng, nx, ny, nz) + + ! Get the index of the diagonals for both rows + call A % search_indices(irow, irow, didx, found) + + ! Perform temporary sums, first do left of diag block, then right of diag block + tmp1 = ZERO + do icol = A % get_row(irow), didx - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) +! print *,A % val(icol),x % val(A % get_col(icol)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) +! print *,A % val(icol),x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) +!print *, irow, b % val(irow), tmp1, A % val(didx), x1 + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 +!stop + end do ROWS + + ! Check convergence + err = sqrt(sum(((tmpx % val - x % val)/tmpx % val)**2)/n) + its = igs + + if (err < tol) exit + + ! Calculation new overrelaxation parameter + w = ONE/(ONE - 0.25_8*cmfd_spectral*w) + + end do GS + + call tmpx % destroy() + + end subroutine cmfd_linsolver_ng + !=============================================================================== ! EXTRACT_RESULTS takes results and puts them in CMFD global data object !=============================================================================== From b93de8ffc6f3f9fa10c54fce4ba5ec8fc6db4836 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 15 Apr 2014 17:29:39 -0400 Subject: [PATCH 11/25] use zero --- src/cmfd_power_solver.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index 65b60bc891..f80c122b71 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -242,7 +242,7 @@ contains subroutine convergence(iter) - use constants, only: ONE, TINY_BIT + use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV @@ -255,7 +255,7 @@ contains kerr = abs(k_o - k_n)/k_n ! Calculate max error in source - where (s_n % val > TINY_BIT) + where (s_n % val > ZERO) serr_v % val = ((s_n % val - s_o % val)/s_n % val)**2 end where serr = sqrt(ONE/dble(s_n % n) * sum(serr_v % val)) From 55d105cf68c881d94121906b0ca21e0eec4d74df Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 18 Apr 2014 10:02:19 -0400 Subject: [PATCH 12/25] fixed using more than 2 groups --- src/cmfd_execute.F90 | 28 ++++++---------------------- 1 file changed, 6 insertions(+), 22 deletions(-) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index fe19bd95b7..645722e618 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -291,7 +291,6 @@ contains logical :: in_mesh ! source site is inside mesh type(StructuredMesh), pointer :: m ! point to mesh - real(8), allocatable :: egrid(:) ! energy grid ! Associate pointer m => meshes(n_user_meshes + 1) @@ -312,19 +311,16 @@ contains cmfd % weightfactors = ONE end if - ! Allocate energy grid and reverse cmfd energy grid - if (.not. allocated(egrid)) allocate(egrid(ng + 1)) - egrid = (/(cmfd % egrid(ng - i + 2), i = 1, ng + 1)/) - ! Compute new weight factors if (new_weights) then ! Zero out weights - cmfd%weightfactors = ZERO + cmfd%weightfactors = ONE - ! Count bank sites in mesh - call count_bank_sites(m, source_bank, cmfd%sourcecounts, egrid, & + ! Count bank sites in mesh and reverse due to egrid structure + call count_bank_sites(m, source_bank, cmfd%sourcecounts, cmfd % egrid, & sites_outside=outside, size_bank=work) + cmfd % sourcecounts = cmfd%sourcecounts(ng:1:-1,:,:,:) ! Check for sites outside of the mesh if (master .and. outside) then @@ -338,25 +334,16 @@ contains cmfd % weightfactors = cmfd % cmfd_src/sum(cmfd % cmfd_src)* & sum(cmfd % sourcecounts) / cmfd % sourcecounts end where -#ifdef CMFD_DEBUG - open(file='source_bank_' // trim(to_str(current_batch)) // '.dat', unit=101) - do i = 1, nx - do j = 1,ny - write(101,*) cmfd % sourcecounts(1,i,j,1) - end do - end do - close(101) -#endif end if -if (.not. cmfd_feedback) return + if (.not. cmfd_feedback) return ! Broadcast weight factors to all procs #ifdef MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & MPI_COMM_WORLD, mpi_err) #endif - end if + end if ! begin loop over source bank do i = 1, int(work,4) @@ -393,9 +380,6 @@ if (.not. cmfd_feedback) return end do - ! Deallocate all - if (allocated(egrid)) deallocate(egrid) - end subroutine cmfd_reweight !=============================================================================== From 47d2fe37ba4598620c1c3134d430d2701e3eed9c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 18 Apr 2014 10:06:34 -0400 Subject: [PATCH 13/25] switch back to 2g solver --- src/cmfd_solver.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index dd343d87fc..8567bd7636 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -162,7 +162,7 @@ contains case(1) cmfd_linsolver => cmfd_linsolver_1g case(2) - cmfd_linsolver => cmfd_linsolver_ng + cmfd_linsolver => cmfd_linsolver_2g case default cmfd_linsolver => cmfd_linsolver_ng end select From b347c46018ddbd8c3903fac2d475ffa4dbd8c82c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 10 Sep 2014 10:00:15 -0400 Subject: [PATCH 14/25] cleaned up code --- src/cmfd_data.F90 | 37 +++++++++++++++++++------------------ src/cmfd_execute.F90 | 3 +-- src/cmfd_input.F90 | 1 + src/cmfd_solver.F90 | 10 ++++++---- src/constants.F90 | 3 --- src/global.F90 | 6 +++--- 6 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index dc5db8352c..ef637b4c0b 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -20,7 +20,8 @@ contains use cmfd_header, only: allocate_cmfd use constants, only: CMFD_NOACCEL - use global, only: cmfd, cmfd_coremap, cmfd_downscatter + use global, only: cmfd, cmfd_coremap, cmfd_downscatter, & + cmfd_fix_balance ! Check for core map and set it up if ((cmfd_coremap) .and. (cmfd%mat_dim == CMFD_NOACCEL)) call set_coremap() @@ -57,12 +58,16 @@ 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 use mesh, only: mesh_indices_to_bin use mesh_header, only: StructuredMesh use string, only: to_str use tally_header, only: TallyObject +#ifdef CMFD_DEBUG + use global, only: current_batch +#endif + 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 @@ -72,7 +77,6 @@ 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 :: ital ! tally object index integer :: ijk(3) ! indices for mesh cell integer :: score_index ! index to pull from tally object @@ -82,10 +86,13 @@ contains 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 +#ifdef CMFD_DEBUG + integer :: l ! iteration counter for leakage debug + real(8) :: leak1 ! group 1 leakage + real(8) :: leak2 ! group 2 leakage +#endif ! Extract spatial and energy indices from object nx = cmfd % indices(1) @@ -709,9 +716,12 @@ contains subroutine compute_dhat() use constants, only: CMFD_NOACCEL, ZERO - use global, only: cmfd, cmfd_coremap, message, dhat_reset, current_batch + use global, only: cmfd, cmfd_coremap, message, dhat_reset use output, only: write_message use string, only: to_str +#ifdef CMFD_DEBUG + use global, only: current_batch +#endif integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -869,8 +879,8 @@ contains function get_reflector_albedo(l, g, i, j, k) - use constants, only: ALBEDO_REJECT - use global, only: cmfd, cmfd_hold_weights + use constants, only: ONE + use global, only: cmfd real(8) :: get_reflector_albedo ! reflector albedo integer, intent(in) :: i ! iteration counter for x @@ -892,8 +902,7 @@ contains ! Calculate albedo if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. & (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then - albedo = ALBEDO_REJECT - cmfd_hold_weights = .true. + albedo = ONE else albedo = (current(2*l-1)/current(2*l))**(shift_idx) end if @@ -958,8 +967,6 @@ contains keff_bal = cmfd % keff_bal - print *, 'KEFF BALANCE is:', keff_bal - ! Return if not two groups if (ng /= 2) return @@ -1007,12 +1014,6 @@ contains nsigf12 = cmfd % nfissxs(1,2,i,j,k) nsigf22 = cmfd % nfissxs(2,2,i,j,k) - ! Check for no fission into group 2 - if (.not.(nsigf12 < 1e-6_8 .and. nsigf22 < 1e-6_8)) then - write(OUTPUT_UNIT,'(A,1PE11.4,1X,1PE11.4)') 'Fission in G=2', & - nsigf12,nsigf22 - end if - ! Compute absorption xs siga1 = sigt1 - sigs11 - sigs12 siga2 = sigt2 - sigs22 - sigs21 diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 645722e618..7aee715ebe 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -84,7 +84,7 @@ contains use global, only: cmfd_begin, cmfd_on, & cmfd_reset, cmfd_run, & - current_batch, cmfd_hold_weights + current_batch ! Check to activate CMFD diffusion and possible feedback ! this guarantees that when cmfd begins at least one batch of tallies are @@ -283,7 +283,6 @@ contains integer :: nz ! maximum number of cells in z direction integer :: ng ! maximum number of energy groups integer :: i ! iteration counter - integer :: j ! iteration counter integer :: ijk(3) ! spatial bin location integer :: e_bin ! energy bin of source particle integer :: n_groups ! number of energy groups diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 0523325c9a..32625034c3 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -66,6 +66,7 @@ contains use xml_interface use, intrinsic :: ISO_FORTRAN_ENV + integer :: i integer :: ng integer, allocatable :: iarray(:) integer, allocatable :: int_array(:) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 8567bd7636..45bc918a14 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -102,9 +102,8 @@ contains use constants, only: ONE, ZERO use error, only: fatal_error - use global, only: cmfd, message, cmfd_shift, keff, & - cmfd_ktol, cmfd_stol - use global, only: cmfd_write_matrices + use global, only: cmfd, cmfd_shift, keff, cmfd_ktol, cmfd_stol, & + cmfd_write_matrices logical :: adjoint @@ -180,10 +179,13 @@ contains subroutine compute_adjoint() use error, only: fatal_error +#ifndef PETSC use global, only: message -#ifdef PETSC +#else use global, only: cmfd_write_matrices +#endif +#ifdef PETSC ! Transpose matrices call loss % transpose() call prod % transpose() diff --git a/src/constants.F90 b/src/constants.F90 index 0d52e8e811..f7d419eaac 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -392,9 +392,6 @@ module constants ! constant to represent a zero flux "albedo" real(8), parameter :: ZERO_FLUX = 999.0_8 - ! constant to represent albedo rejection - real(8), parameter :: ALBEDO_REJECT = 999.0_8 - ! constant for writing out no residual real(8), parameter :: CMFD_NORES = 99999.0_8 diff --git a/src/global.F90 b/src/global.F90 index 3b9b8bfeea..e2509ce1a7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -325,6 +325,9 @@ module global ! Flag to reset dhats to zero logical :: dhat_reset = .false. + ! Flag to enforce neutron balance + logical :: cmfd_fix_balance = .false. + ! Flag to activate neutronic feedback via source weights logical :: cmfd_feedback = .false. @@ -332,9 +335,6 @@ module global integer :: n_cmfd_meshes = 1 ! # of structured meshes integer :: n_cmfd_tallies = 3 ! # of user-defined tallies - ! Flag to hold cmfd weight adjustment factors - logical :: cmfd_hold_weights = .false. - ! Eigenvalue solver type character(len=10) :: cmfd_solver_type = 'power' From bb1785c0629dd65c5fca7c439f3430d7f447a42e Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 16:51:12 -0400 Subject: [PATCH 15/25] removed subroutine to enforce neutron balance --- src/cmfd_data.F90 | 148 +-------------------------------------------- src/cmfd_input.F90 | 8 --- src/global.F90 | 3 - 3 files changed, 1 insertion(+), 158 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ef637b4c0b..116c9bacb1 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -20,8 +20,7 @@ contains use cmfd_header, only: allocate_cmfd use constants, only: CMFD_NOACCEL - use global, only: cmfd, cmfd_coremap, cmfd_downscatter, & - cmfd_fix_balance + use global, only: cmfd, cmfd_coremap, cmfd_downscatter ! Check for core map and set it up if ((cmfd_coremap) .and. (cmfd%mat_dim == CMFD_NOACCEL)) call set_coremap() @@ -32,9 +31,6 @@ contains ! Compute effective downscatter cross section if (cmfd_downscatter) call compute_effective_downscatter() - ! Compute perfect balance - if (cmfd_fix_balance) call fix_neutron_balance() - ! Check neutron balance call neutron_balance() @@ -912,148 +908,6 @@ contains end function get_reflector_albedo -!=============================================================================== -! FIX_NEUTRON_BALANCE is a method to adjust parameters to have perfect balance -!=============================================================================== - - subroutine fix_neutron_balance() - - use constants, only: ONE, ZERO, CMFD_NOACCEL - use global, only: cmfd, message - use output, only: write_message - - 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 :: l ! iteration counter for surface - real(8) :: leak1 ! leakage rate in group 1 - real(8) :: leak2 ! leakage rate in group 2 - real(8) :: flux1 ! group 1 volume int flux - real(8) :: flux2 ! group 2 volume int flux - real(8) :: sigt1 ! group 1 total xs - real(8) :: sigt2 ! group 2 total xs - real(8) :: sigs11 ! scattering transfer 1 --> 1 - real(8) :: sigs21 ! scattering transfer 2 --> 1 - real(8) :: sigs12 ! scattering transfer 1 --> 2 - real(8) :: sigs22 ! scattering transfer 2 --> 2 - real(8) :: nsigf11 ! fission transfer 1 --> 1 - real(8) :: nsigf21 ! fission transfer 2 --> 1 - real(8) :: nsigf12 ! fission transfer 1 --> 2 - real(8) :: nsigf22 ! fission transfer 2 --> 2 - real(8) :: siga1 ! group 1 abs xs - real(8) :: siga2 ! group 2 abs xs - real(8) :: sigs12_eff ! effective downscatter xs - real(8) :: a ! matrix (1,1) element for balance - real(8) :: b ! matrix (1,2) element for balance - real(8) :: c ! matrix (2,1) element for balance - real(8) :: d ! matrix (2,2) element for balance - real(8) :: r1 ! right hand side element 1 - real(8) :: r2 ! right hand side element 2 - real(8) :: det ! determinant of balance matrix - real(8) :: keff_bal ! keffective for balance eq. - - message = 'Correcting neutron balance' - call write_message(1) - - ! Extract spatial and energy indices from object - nx = cmfd % indices(1) - ny = cmfd % indices(2) - nz = cmfd % indices(3) - ng = cmfd % indices(4) - - keff_bal = cmfd % keff_bal - - ! Return if not two groups - if (ng /= 2) return - - ! Begin loop around space and energy groups - ZLOOP: do k = 1, nz - - YLOOP: do j = 1, ny - - XLOOP: do i = 1, nx - - ! Check for active mesh - if (allocated(cmfd%coremap)) then - if (cmfd%coremap(i,j,k) == CMFD_NOACCEL) cycle - end if - - ! Compute leakage in groups 1 and 2 - leak1 = ZERO - leak2 = ZERO - LEAK: do l = 1, 3 - - leak1 = leak1 + ((cmfd % current(4*l,1,i,j,k) - & - cmfd % current(4*l-1,1,i,j,k))) - & - ((cmfd % current(4*l-2,1,i,j,k) - & - cmfd % current(4*l-3,1,i,j,k))) - - leak2 = leak2 + ((cmfd % current(4*l,2,i,j,k) - & - cmfd % current(4*l-1,2,i,j,k))) - & - ((cmfd % current(4*l-2,2,i,j,k) - & - cmfd % current(4*l-3,2,i,j,k))) - - - end do LEAK - - ! Extract cross sections and flux from object - flux1 = cmfd % flux(1,i,j,k) - flux2 = cmfd % flux(2,i,j,k) - sigt1 = cmfd % totalxs(1,i,j,k) - sigt2 = cmfd % totalxs(2,i,j,k) - sigs11 = cmfd % scattxs(1,1,i,j,k) - sigs21 = cmfd % scattxs(2,1,i,j,k) - sigs12 = cmfd % scattxs(1,2,i,j,k) - sigs22 = cmfd % scattxs(2,2,i,j,k) - nsigf11 = cmfd % nfissxs(1,1,i,j,k) - nsigf21 = cmfd % nfissxs(2,1,i,j,k) - nsigf12 = cmfd % nfissxs(1,2,i,j,k) - nsigf22 = cmfd % nfissxs(2,2,i,j,k) - - ! Compute absorption xs - siga1 = sigt1 - sigs11 - sigs12 - siga2 = sigt2 - sigs22 - sigs21 - - ! Create matrix and solve for effective downscatter and thermal flux - a = flux1 - b = -ONE/keff_bal*nsigf21 - c = flux1 - d = -siga2 + ONE/keff_bal*nsigf22 - det = a*d - b*c - r1 = ONE/keff_bal*nsigf11*flux1 - siga1*flux1 - leak1 - r2 = leak2 - ONE/keff_bal*nsigf12*flux1 - sigs12_eff = ONE/det*(d*r1 - b*r2) - flux2 = ONE/det*(a*r2 - c*r1) - - ! Recompute total cross sections (use effective and no upscattering) - sigt1 = siga1 + sigs11 + sigs12_eff - sigt2 = siga2 + sigs22 - - ! Record total xs - cmfd % totalxs(1,i,j,k) = sigt1 - cmfd % totalxs(2,i,j,k) = sigt2 - - ! Record effective downscatter xs - cmfd % scattxs(1,2,i,j,k) = sigs12_eff - - ! Zero out upscatter cross section - cmfd % scattxs(2,1,i,j,k) = ZERO - - ! Record thermal flux - cmfd % flux(2,i,j,k) = flux2 - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine fix_neutron_balance - !=============================================================================== ! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter !=============================================================================== diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 32625034c3..2e531a8157 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -174,14 +174,6 @@ contains dhat_reset = .true. end if - ! Enforce perfect neutron balance - if (check_for_node(doc, "balance")) then - call get_node_value(doc, "balance", temp_str) - call lower_case(temp_str) - if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_fix_balance = .true. - end if - ! Set the solver type if (check_for_node(doc, "solver")) & call get_node_value(doc, "solver", cmfd_solver_type) diff --git a/src/global.F90 b/src/global.F90 index e2509ce1a7..cad98dc329 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -325,9 +325,6 @@ module global ! Flag to reset dhats to zero logical :: dhat_reset = .false. - ! Flag to enforce neutron balance - logical :: cmfd_fix_balance = .false. - ! Flag to activate neutronic feedback via source weights logical :: cmfd_feedback = .false. From f1f7e76d60fa36137d01e560cfb0103ed3f31e10 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 16:55:50 -0400 Subject: [PATCH 16/25] removed CMFD debugging code --- src/CMakeLists.txt | 10 ----- src/cmfd_data.F90 | 101 ------------------------------------------- src/cmfd_execute.F90 | 9 ---- 3 files changed, 120 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3f336767cd..8a9424e79d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -26,7 +26,6 @@ option(debug "Compile with debug flags" OFF) option(optimize "Turn on all compiler optimization flags" OFF) option(verbose "Create verbose Makefiles" OFF) option(coverage "Compile with flags" OFF) -option(cmfd_debug "Compile with cmfd debug flags" OFF) if (verbose) set(CMAKE_VERBOSE_MAKEFILE on) @@ -215,15 +214,6 @@ if(GIT_SHA1_SUCCESS EQUAL 0) add_definitions(-DGIT_SHA1="${GIT_SHA1}") endif() -#=============================================================================== -# CMFD debug flags -#=============================================================================== - -if (cmfd_debug) - message("-- CMFD Debug flags on") - add_definitions(-DCMFD_DEBUG) -endif() - #=============================================================================== # FoX Fortran XML Library #=============================================================================== diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 116c9bacb1..3ed77d78be 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -60,10 +60,6 @@ contains use string, only: to_str use tally_header, only: TallyObject -#ifdef CMFD_DEBUG - use global, only: current_batch -#endif - 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 @@ -84,11 +80,6 @@ contains real(8) :: flux ! temp variable for flux type(TallyObject), pointer :: t => null() ! pointer for tally object type(StructuredMesh), pointer :: m => null() ! pointer for mesh object -#ifdef CMFD_DEBUG - integer :: l ! iteration counter for leakage debug - real(8) :: leak1 ! group 1 leakage - real(8) :: leak2 ! group 2 leakage -#endif ! Extract spatial and energy indices from object nx = cmfd % indices(1) @@ -311,85 +302,6 @@ contains if (associated(t)) nullify(t) if (associated(m)) nullify(m) -#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) - open(file='totalxs2_' // trim(to_str(current_batch)) // '.dat', unit=104) - open(file='p1scattxs1_' // trim(to_str(current_batch)) // '.dat', unit=105) - open(file='p1scattxs2_' // trim(to_str(current_batch)) // '.dat', unit=106) - open(file='scattxs11_' // trim(to_str(current_batch)) // '.dat', unit=107) - open(file='scattxs12_' // trim(to_str(current_batch)) // '.dat', unit=108) - open(file='scattxs21_' // trim(to_str(current_batch)) // '.dat', unit=109) - open(file='scattxs22_' // trim(to_str(current_batch)) // '.dat', unit=110) - open(file='nufissxs11_' // trim(to_str(current_batch)) // '.dat', unit=111) - open(file='nufissxs12_' // trim(to_str(current_batch)) // '.dat', unit=112) - open(file='nufissxs21_' // trim(to_str(current_batch)) // '.dat', unit=113) - open(file='nufissxs22_' // trim(to_str(current_batch)) // '.dat', unit=114) - open(file='diff_coef1_' // trim(to_str(current_batch)) // '.dat', unit=115) - open(file='diff_coef2_' // trim(to_str(current_batch)) // '.dat', unit=116) - open(file='flux1_' // trim(to_str(current_batch)) // '.dat', unit=117) - open(file='flux2_' // trim(to_str(current_batch)) // '.dat', unit=118) - open(file='leak1_' // trim(to_str(current_batch)) // '.dat', unit=119) - open(file='leak2_' // trim(to_str(current_batch)) // '.dat', unit=120) - - do i = 1, nx - do j = 1,ny - leak1 = ZERO - leak2 = ZERO - do l = 1, 3 - leak1 = leak1 + ((cmfd % current(4*l,1,i,j,1) - & - cmfd % current(4*l-1,1,i,j,1))) - & - ((cmfd % current(4*l-2,1,i,j,1) - & - cmfd % current(4*l-3,1,i,j,1))) - end do - do l = 1, 3 - leak2 = leak2 + ((cmfd % current(4*l,2,i,j,1) - & - cmfd % current(4*l-1,2,i,j,1))) - & - ((cmfd % current(4*l-2,2,i,j,1) - & - cmfd % current(4*l-3,2,i,j,1))) - end do - write(102,*) cmfd % openmc_src(1,i,j,1) - write(103,*) cmfd % totalxs(1,i,j,1) - write(104,*) cmfd % totalxs(2,i,j,1) - write(105,*) cmfd % p1scattxs(1,i,j,1) - write(106,*) cmfd % p1scattxs(2,i,j,1) - write(107,*) cmfd % scattxs(1,1,i,j,1) - write(108,*) cmfd % scattxs(1,2,i,j,1) - write(109,*) cmfd % scattxs(2,1,i,j,1) - write(110,*) cmfd % scattxs(2,2,i,j,1) - write(111,*) cmfd % nfissxs(1,1,i,j,1) - write(112,*) cmfd % nfissxs(1,2,i,j,1) - write(113,*) cmfd % nfissxs(2,1,i,j,1) - write(114,*) cmfd % nfissxs(2,2,i,j,1) - write(115,*) cmfd % diffcof(1,i,j,1) - write(116,*) cmfd % diffcof(2,i,j,1) - write(117,*) cmfd % flux(1,i,j,1) - write(118,*) cmfd % flux(2,i,j,1) - write(119,*) leak1 - write(120,*) leak2 - end do - end do - close(102) - close(103) - close(104) - close(105) - close(106) - close(107) - close(108) - close(109) - close(110) - close(111) - close(112) - close(113) - close(114) - close(115) - close(116) - close(117) - close(118) - close(119) - close(120) -#endif - end subroutine compute_xs !=============================================================================== @@ -715,9 +627,6 @@ contains use global, only: cmfd, cmfd_coremap, message, dhat_reset use output, only: write_message use string, only: to_str -#ifdef CMFD_DEBUG - use global, only: current_batch -#endif integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -751,9 +660,6 @@ contains nxyz(1,:) = (/1,nx/) nxyz(2,:) = (/1,ny/) nxyz(3,:) = (/1,nz/) -#ifdef CMFD_DEBUG - open(file='cmfd_dhat_' // trim(to_str(current_batch)) // '.dat', unit=125) -#endif ! Geting loop over group and spatial indices ZLOOP: do k = 1,nz @@ -844,9 +750,6 @@ contains cmfd%dhat(l,g,i,j,k) = ZERO end if -#ifdef CMFD_DEBUG - write(125,*) dhat -#endif end do LEAK end do GROUP @@ -863,10 +766,6 @@ contains call write_message(1) end if -#ifdef CMFD_DEBUG - close(unit=125) -#endif - end subroutine compute_dhat !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 7aee715ebe..f76dfb42f9 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -238,15 +238,6 @@ contains cmfd % src_cmp(current_batch) = sqrt(ONE/cmfd % norm * & sum((cmfd % cmfd_src - cmfd % openmc_src)**2)) -#ifdef CMFD_DEBUG - open(file='cmfd_src_' // trim(to_str(current_batch)) // '.dat', unit=100) - do i = 1, nx - do j = 1,ny - write(100,*) cmfd % cmfd_src(1,i,j,1) - end do - end do - close(100) -#endif end if #ifdef MPI From 0953dc5edf6ae516583f5b212d53d6d773bf572c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 16:56:40 -0400 Subject: [PATCH 17/25] removed balance from users guide --- docs/source/usersguide/input.rst | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 1a8d1f2ac7..5a775be8e9 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1312,18 +1312,6 @@ and not when PETSc is active. *Default*: 1.e-10 -```` Element ---------------------- - -The ```` element controls whether exact neutron balance should be enforced -from CMFD tallies before creating CMFD matrices. This changes effective downscatter -cross section and thermal flux to create exact balance. It should be noted that -this option has led to instabilities when performing CMFD. It can be turned on -with "true" and off with "false". - - *Default*: false - - ```` Element ------------------- From 12f9f077573e6177eb0790d06f731c83c8bc995a Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 16:57:37 -0400 Subject: [PATCH 18/25] modified comment on default weightfactors value --- src/cmfd_execute.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index f76dfb42f9..d14a3af243 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -304,7 +304,7 @@ contains ! Compute new weight factors if (new_weights) then - ! Zero out weights + ! Set weight factors to a default 1.0 cmfd%weightfactors = ONE ! Count bank sites in mesh and reverse due to egrid structure From 495bab1209e0d16f19b1db8845301e84102d8156 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 17:09:05 -0400 Subject: [PATCH 19/25] fixed compilation error, added check for max GS iters --- src/cmfd_data.F90 | 4 +++- src/cmfd_input.F90 | 2 +- src/cmfd_solver.F90 | 43 +++++++++++++++++++++++++++++++------------ src/input_xml.F90 | 2 +- 4 files changed, 36 insertions(+), 15 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 3ed77d78be..abba2d47d7 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -211,7 +211,9 @@ contains ! Bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & t % results(2,score_index) % sum - cmfd % keff_bal = cmfd % keff_bal + t % results(2,score_index) % sum / dble(t % n_realizations) + cmfd % keff_bal = cmfd % keff_bal + & + t % results(2,score_index) % sum / & + dble(t % n_realizations) end do INGROUP diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 98aaed9c7d..22bb7b1c8f 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -169,7 +169,7 @@ contains ! Reset dhat parameters if (check_for_node(doc, "dhat_reset")) then call get_node_value(doc, "dhat_reset", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & dhat_reset = .true. end if diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 38e3436a24..df19861e57 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -356,7 +356,8 @@ contains subroutine cmfd_linsolver_1g(A, b, x, tol, its) use constants, only: ONE, ZERO - use global, only: cmfd, cmfd_spectral + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message type(Matrix), intent(inout) :: A ! coefficient matrix type(Vector), intent(inout) :: b ! right hand side vector @@ -398,6 +399,12 @@ contains ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + ! Copy over x vector call tmpx % copy(x) @@ -456,7 +463,8 @@ contains subroutine cmfd_linsolver_2g(A, b, x, tol, its) use constants, only: ONE, ZERO - use global, only: cmfd, cmfd_spectral + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message type(Matrix), intent(inout) :: A ! coefficient matrix type(Vector), intent(inout) :: b ! right hand side vector @@ -510,6 +518,12 @@ contains ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + ! Copy over x vector call tmpx % copy(x) @@ -595,13 +609,14 @@ contains subroutine cmfd_linsolver_ng(A, b, x, tol, its) use constants, only: ONE, ZERO - use global, only: cmfd, cmfd_spectral + use error, only: fatal_error + use global, only: cmfd, cmfd_spectral, message - type(Matrix) :: A ! coefficient matrix - type(Vector) :: b ! right hand side vector - type(Vector) :: x ! unknown vector - real(8) :: tol ! tolerance on final error - integer :: its ! number of inner iterations + type(Matrix), intent(inout) :: A ! coefficient matrix + type(Vector), intent(inout) :: b ! right hand side vector + type(Vector), intent(inout) :: x ! unknown vector + real(8), intent(in) :: tol ! tolerance on final error + integer, intent(out) :: its ! number of inner iterations integer :: g ! group index integer :: i ! loop counter for x @@ -636,6 +651,12 @@ contains ! Perform Gauss Seidel iterations GS: do igs = 1, 10000 + ! Check for max iterations met + if (igs == 10000) then + message = 'Maximum Gauss-Seidel iterations encountered.' + call fatal_error() + endif + ! Copy over x vector call tmpx % copy(x) @@ -652,19 +673,17 @@ contains tmp1 = ZERO do icol = A % get_row(irow), didx - 1 tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) -! print *,A % val(icol),x % val(A % get_col(icol)) end do do icol = didx + 1, A % get_row(irow + 1) - 1 tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) -! print *,A % val(icol),x % val(A % get_col(icol)) end do ! Solve for new x x1 = (b % val(irow) - tmp1)/A % val(didx) -!print *, irow, b % val(irow), tmp1, A % val(didx), x1 + ! Perform overrelaxation x % val(irow) = (ONE - w)*x % val(irow) + w*x1 -!stop + end do ROWS ! Check convergence diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5d0a159238..c7bf3ef17a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -268,7 +268,7 @@ contains ! Check if we want to write out source if (check_for_node(node_source, "write_initial")) then call get_node_value(node_source, "write_initial", temp_str) - call lower_case(temp_str) + temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & write_initial_source = .true. end if From 3b1d04eeaf0a53ca4a4dc926169d3be1594eb4bd Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 17:14:54 -0400 Subject: [PATCH 20/25] added input gauss_seidel_tolerance to replace atoli and rtoli --- src/cmfd_input.F90 | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 22bb7b1c8f..235a23c60b 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -68,12 +68,14 @@ contains integer :: i integer :: ng + integer :: n_params integer, allocatable :: iarray(:) integer, allocatable :: int_array(:) logical :: file_exists ! does cmfd.xml exist? logical :: found character(MAX_LINE_LEN) :: filename character(MAX_LINE_LEN) :: temp_str + real(8) :: gs_tol(2) type(Node), pointer :: doc => null() type(Node), pointer :: node_mesh => null() @@ -256,10 +258,17 @@ contains call get_node_value(doc, "ktol", cmfd_ktol) if (check_for_node(doc, "stol")) & call get_node_value(doc, "stol", cmfd_stol) - if (check_for_node(doc, "atoli")) & - call get_node_value(doc, "atoli", cmfd_atoli) - if (check_for_node(doc, "rtoli")) & - call get_node_value(doc, "rtoli", cmfd_rtoli) + if (check_for_node(doc, "gauss_seidel_tolerance")) then + n_params = get_arraysize_integer(doc, "gauss_seidel_tolerance") + if (n_params /= 2) then + message = 'Gauss Seidel tolerance is not 2 parameters & + &(absolute, relative).' + call fatal_error() + end if + call get_node_array(doc, "gauss_seidel_tolerance", gs_tol) + cmfd_atoli = gs_tol(1) + cmfd_rtoli = gs_tol(2) + end if ! Create tally objects call create_cmfd_tally(doc) From 364199eb2847d67e15d8bc3544a9038b5b3ce177 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 17:25:06 -0400 Subject: [PATCH 21/25] added input description of write_initial_source and allowed binary extension for name --- docs/source/usersguide/input.rst | 42 ++++++++++++++++---------------- src/source.F90 | 4 +++ 2 files changed, 25 insertions(+), 21 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 5a775be8e9..4fc49b35f5 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -433,6 +433,13 @@ attributes/sub-elements: *Default*: 0.988 2.249 + :write_initial: + An element specifying whether to write out the intial source bank used + at the beginning of the first batch. The output file is named + "initial_source.binary(h5)" + + *Default*: false + ```` Element ------------------------- @@ -1303,15 +1310,6 @@ Currently, it allows users to accelerate fission source convergence during inactive neutron batches. To run CMFD, the ```` element in ``settings.xml`` should be set to "true". -```` Element --------------------- - -The ```` element specifies the absolute inner tolerance on the Gauss-Seidel -when performing CMFD calculations. It is only used in the standalone CMFD solver -and not when PETSc is active. - - *Default*: 1.e-10 - ```` Element ------------------- @@ -1338,9 +1336,9 @@ The ```` element sets one additional CMFD output column. Options are: ```` Element ------------------------ -The ```` element controls whether dhat nonlinear CMFD parameters -should be reset to zero before solving CMFD eigenproblem. It can be turned on with -"true" and off with "false". +The ```` element controls whether :math:`\widehat{D}` nonlinear +CMFD parameters should be reset to zero before solving CMFD eigenproblem. +It can be turned on with "true" and off with "false". *Default*: false @@ -1362,6 +1360,17 @@ It can be turned on with "true" and off with "false". *Default*: false +```` Element +-------------------- + +The ```` element specifies two parameters. The first is +the absolute inner tolerance for Gauss-Seidel iterations when performing CMFD +and the second is the relative inner tolerance for Gauss-Seidel iterations +for CMFD calculations. It is only used in the standalone CMFD power iteration +solver and not when PETSc is active. + + *Default*: 1.e-10 1.e-5 + ```` Element ------------------------- @@ -1455,15 +1464,6 @@ This option can be turned on with "true" and turned off with "false". *Default*: false -```` Element --------------------- - -The ```` element specifies the relative inner tolerance on the Gauss-Seidel -when performing CMFD calculations. It is only used in the standalone CMFD solver -and not when PETSc is active. - - *Default*: 1.e-5 - ```` Element ------------------------- diff --git a/src/source.F90 b/src/source.F90 index f6de6042b1..fd889abc42 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -81,7 +81,11 @@ contains if (write_initial_source) then message = 'Writing out initial source guess...' call write_message(1) +#ifdef HDF5 filename = 'initial_source.h5' +#else + filename = 'initial_source.binary' +#endif call sp % file_create(filename, serial = .false.) call sp % write_source_bank() call sp % file_close() From 5dab31aa49830907ba981ba85e8cc6ca901117b2 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 17:26:31 -0400 Subject: [PATCH 22/25] put path_output in front of initial source file --- src/source.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/source.F90 b/src/source.F90 index fd889abc42..9575020e6d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -82,9 +82,9 @@ contains message = 'Writing out initial source guess...' call write_message(1) #ifdef HDF5 - filename = 'initial_source.h5' + filename = trim(path_output) // 'initial_source.h5' #else - filename = 'initial_source.binary' + filename = trim(path_output) // 'initial_source.binary' #endif call sp % file_create(filename, serial = .false.) call sp % write_source_bank() From a531aea8ef9b448b89fdd0b83827ca502f1e6e72 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 19:24:00 -0400 Subject: [PATCH 23/25] changed test input files to use new gs_tol input --- tests/test_cmfd_feed/cmfd.xml | 4 +--- tests/test_cmfd_nofeed/cmfd.xml | 3 +-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/test_cmfd_feed/cmfd.xml b/tests/test_cmfd_feed/cmfd.xml index 390e558639..228b260a46 100644 --- a/tests/test_cmfd_feed/cmfd.xml +++ b/tests/test_cmfd_feed/cmfd.xml @@ -12,7 +12,5 @@ dominance power true - 1.e-15 - 1.e-20 - + 1.e-15 1.e-20 diff --git a/tests/test_cmfd_nofeed/cmfd.xml b/tests/test_cmfd_nofeed/cmfd.xml index 4a4bf30108..eb6c2c721a 100644 --- a/tests/test_cmfd_nofeed/cmfd.xml +++ b/tests/test_cmfd_nofeed/cmfd.xml @@ -12,7 +12,6 @@ dominance power false - 1.e-15 - 1.e-20 + 1.e-15 1.e-20 From 87d7a640b02939764eac2453bea001a2e04f2b73 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 19:45:40 -0400 Subject: [PATCH 24/25] gs tol should be read in as double --- src/cmfd_input.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 235a23c60b..67e8f0d0f3 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -259,7 +259,7 @@ contains if (check_for_node(doc, "stol")) & call get_node_value(doc, "stol", cmfd_stol) if (check_for_node(doc, "gauss_seidel_tolerance")) then - n_params = get_arraysize_integer(doc, "gauss_seidel_tolerance") + n_params = get_arraysize_double(doc, "gauss_seidel_tolerance") if (n_params /= 2) then message = 'Gauss Seidel tolerance is not 2 parameters & &(absolute, relative).' From 855483e540ecf64308d8418abf3ac37715c394f3 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 17 Sep 2014 19:56:18 -0400 Subject: [PATCH 25/25] updated RELAX NG schemata --- src/relaxng/cmfd.rnc | 25 +++++++++++++------------ src/relaxng/settings.rnc | 3 ++- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/relaxng/cmfd.rnc b/src/relaxng/cmfd.rnc index cbcb7709ae..bbb7fc9e6b 100644 --- a/src/relaxng/cmfd.rnc +++ b/src/relaxng/cmfd.rnc @@ -22,15 +22,9 @@ element cmfd { element feedback { xsd:boolean }? & - element n_cmfd_procs { xsd:int }? & - - element reset { xsd:boolean }? & - - element balance { xsd:boolean }? & - element downscatter { xsd:boolean }? & - element run_2grp { xsd:boolean }? & + element dhat_reset { xsd:boolean }? & element solver { xsd:string }? & @@ -40,8 +34,6 @@ element cmfd { element power_monitor { xsd:boolean }? & - element write_balance { xsd:boolean }? & - element write_matrices { xsd:boolean }? & element run_adjoint { xsd:boolean }? & @@ -50,9 +42,18 @@ element cmfd { element begin { xsd:int }? & - element inactive { xsd:boolean }? & + element tally_reset { list { xsd:int+ } }? & - element active_flush { xsd:int }? & + element display { xsd:string }? & + + element spectral { xsd:double }? & + + element shift { xsd:double }? & + + element ktol { xsd: double }? & + + element stol { xsd: double }? & + + element gauss_seidel_tolerance { list { xsd:double+ } }? - element keff_tol { xsd:double }? } diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index dd4b91dc7d..707ecbbc80 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -85,7 +85,8 @@ element settings { attribute interplation { xsd:string { maxLength = "10" } })? & (element parameters { list { xsd:double+ } } | attribute parameters { list { xsd:double+ } })? - }? + }? & + (element write_initial { xsd:boolean } | attribute write_initial { xsd:boolean })? }? & element state_point {