diff --git a/docs/source/releasenotes/index.rst b/docs/source/releasenotes/index.rst index b40126b243..9799a5cfc7 100644 --- a/docs/source/releasenotes/index.rst +++ b/docs/source/releasenotes/index.rst @@ -10,6 +10,7 @@ bugs fixed, and known issues for each successive release. .. toctree:: :maxdepth: 1 + notes_0.6.1 notes_0.6.0 notes_0.5.4 notes_0.5.3 diff --git a/docs/source/releasenotes/notes_0.6.1.rst b/docs/source/releasenotes/notes_0.6.1.rst new file mode 100644 index 0000000000..b0651a8629 --- /dev/null +++ b/docs/source/releasenotes/notes_0.6.1.rst @@ -0,0 +1,63 @@ +.. _notes_0.6.1: + +============================== +Release Notes for OpenMC 0.6.1 +============================== + +------------------- +System Requirements +------------------- + +There are no special requirements for running the OpenMC code. As of this +release, OpenMC has been tested on a variety of Linux distributions, Mac OS X, +and Microsoft Windows 7. Memory requirements will vary depending on the size of +the problem at hand (mostly on the number of nuclides in the problem). + +------------ +New Features +------------ + +- Coarse mesh finite difference acceleration no longer requires PETSc +- Statepoint file numbering is now zero-padded +- Python scripts now compatible with Python 2 or 3 +- Ability to run particle restarts in fixed source calculations +- Capability to filter box source by fissionable materials +- Nuclide/element names are now case insensitive in input files +- Improved treatment of resonance scattering for heavy nuclides + +--------- +Bug Fixes +--------- + +- 03e890_: Check for energy-dependent multiplicities in ACE files +- 4439de_: Fix distance-to-surface calculation for general plane surface +- 5808ed_: Account for differences in URR band probabilities at different energies +- 2e60c0_: Allow zero atom/weight percents in materials +- 3e0870_: Don't use PWD environment variable when setting path to input files +- dc4776_: Handle probability table resampling correctly +- 01178b_: Fix metastables nuclides in NNDC cross_sections.xml file +- 62ec43_: Don't read tallies.xml when OpenMC is run in plotting mode +- 2a95ef_: Prevent segmentation fault on "current" score without mesh filter + +.. _03e890: https://github.com/mit-crpg/openmc/commit/03e890 +.. _4439de: https://github.com/mit-crpg/openmc/commit/4439de +.. _5808ed: https://github.com/mit-crpg/openmc/commit/5808ed +.. _2e60c0: https://github.com/mit-crpg/openmc/commit/2e60c0 +.. _3e0870: https://github.com/mit-crpg/openmc/commit/3e0870 +.. _dc4776: https://github.com/mit-crpg/openmc/commit/dc4776 +.. _01178b: https://github.com/mit-crpg/openmc/commit/01178b +.. _62ec43: https://github.com/mit-crpg/openmc/commit/62ec43 +.. _2a95ef: https://github.com/mit-crpg/openmc/commit/2a95ef + +------------ +Contributors +------------ + +This release contains new contributions from the following people: + +- `Sterling Harper `_ +- `Bryan Herman `_ +- `Adam Nelson `_ +- `Paul Romano `_ +- `Jon Walsh `_ +- `Will Boyd `_ diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index de3cbc7f23..cbe06e5a0d 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -261,7 +261,7 @@ or sub-elements and can be set to either "false" or "true". *Default*: true ```` Element ----------------------- +---------------------------------- The ``resonance_scattering`` element can contain one or more of the following attributes or sub-elements: @@ -269,7 +269,7 @@ attributes or sub-elements: :scatterer: An element with attributes/sub-elements called ``nuclide``, ``method``, ``xs_label``, ``xs_label_0K``, ``E_min``, and ``E_max``. The ``nuclide`` - attribute is the name, as given by the ``name`` attribute within the + attribute is the name, as given by the ``name`` attribute within the ``nuclide`` sub-element of the ``material`` element in ``materials.xml``, of the nuclide to which a resonance scattering treatment is to be applied. The ``method`` attribute gives the type of resonance scattering treatment @@ -433,6 +433,13 @@ attributes/sub-elements: *Default*: 0.988 2.249 + :write_initial: + An element specifying whether to write out the initial 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 controls the batch where CMFD tallies should be -reset. CMFD tallies should be reset before active batches so they are accumulated -without bias. - - *Default*: 0 - ```` Element ------------------- @@ -1333,7 +1331,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 :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 + +```` 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 ---------------------- @@ -1344,24 +1360,16 @@ It can be turned on with "true" and off with "false". *Default*: false -```` Element ----------------------- +```` 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" +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*: 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 + *Default*: 1.e-10 1.e-5 ```` Element ------------------------- @@ -1370,9 +1378,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 ------------------ @@ -1441,14 +1456,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 --------------------------- @@ -1461,16 +1468,8 @@ This option can be turned on with "true" and turned off with "false". ------------------------- 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 @@ -1483,6 +1482,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_data.F90 b/src/cmfd_data.F90 index b5612dcd05..abba2d47d7 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 !============================================================================== @@ -103,6 +101,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 @@ -211,6 +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) end do INGROUP @@ -623,7 +626,9 @@ 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 + use string, only: to_str integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -743,7 +748,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 @@ -755,6 +762,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 !=============================================================================== @@ -763,8 +776,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 @@ -786,8 +799,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 @@ -797,137 +809,6 @@ contains end function get_reflector_albedo -!=============================================================================== -! 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 - - 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 - - ! Extract spatial and energy indices from object - nx = cmfd % indices(1) - ny = cmfd % indices(2) - nz = cmfd % indices(3) - ng = cmfd % indices(4) - - ! 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) - - ! 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 - - ! 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 - - ! 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 - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine fix_neutron_balance -#endif - !=============================================================================== ! COMPUTE_EFFECTIVE_DOWNSCATTER changes downscatter rate for zero upscatter !=============================================================================== diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 35485a8945..d14a3af243 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() @@ -82,36 +82,23 @@ contains subroutine cmfd_init_batch() - use global, only: cmfd_begin, cmfd_on, cmfd_tally_on, & - cmfd_inact_flush, cmfd_act_flush, cmfd_run, & - current_batch, cmfd_hold_weights + use global, only: cmfd_begin, cmfd_on, & + cmfd_reset, cmfd_run, & + current_batch ! Check to activate CMFD diffusion and possible feedback ! this guarantees that when cmfd begins at least one batch of tallies are ! 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 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 @@ -144,6 +131,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 @@ -272,6 +260,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 @@ -292,7 +281,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) @@ -313,19 +301,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 + ! Set weight factors to a default 1.0 + 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 @@ -341,12 +326,14 @@ contains end where 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, & MPI_COMM_WORLD, mpi_err) #endif - end if + end if ! begin loop over source bank do i = 1, int(work,4) @@ -383,9 +370,6 @@ contains end do - ! Deallocate all - if (allocated(egrid)) deallocate(egrid) - end subroutine cmfd_reweight !=============================================================================== 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 diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 5ab4caf480..67e8f0d0f3 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -66,12 +66,16 @@ contains use xml_interface use, intrinsic :: ISO_FORTRAN_ENV + 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() @@ -164,6 +168,14 @@ contains cmfd_downscatter = .true. end if + ! Reset dhat parameters + if (check_for_node(doc, "dhat_reset")) then + call get_node_value(doc, "dhat_reset", temp_str) + temp_str = to_lower(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) @@ -201,30 +213,31 @@ contains call get_node_value(doc, "run_adjoint", temp_str) temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_run_adjoint = .true. +#ifndef PETSC + message = 'Must use PETSc when running adjoint option.' + call fatal_error() +#endif + cmfd_run_adjoint = .true. end if ! Batch to begin cmfd 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) - temp_str = to_lower(temp_str) - if (trim(temp_str) == 'false' .or. trim(temp_str) == '0') & - cmfd_tally_on = .false. + ! 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 - - ! 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) ! Get display if (check_for_node(doc, "display")) & @@ -245,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_double(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) 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)) diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 69a2829039..b3e8fb8a88 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -9,7 +9,7 @@ module cmfd_solver implicit none private - public :: cmfd_solver_execute + public :: cmfd_solver_execute real(8) :: k_n ! new k-eigenvalue real(8) :: k_o ! old k-eigenvalue @@ -42,7 +42,7 @@ module cmfd_solver type(Vector), intent(inout) :: x real(8), intent(in) :: tol integer, intent(out) :: i - end subroutine linsolve + end subroutine linsolve end interface contains @@ -80,7 +80,7 @@ contains ! Stop timer for build call time_cmfdbuild % stop() - ! Begin power iteration + ! Begin power iteration call time_cmfdsolve % start() call execute_power_iter() call time_cmfdsolve % stop() @@ -88,7 +88,7 @@ contains ! Extract results call extract_results() - ! Deallocate data + ! Deallocate data call finalize() end subroutine cmfd_solver_execute @@ -101,9 +101,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, intent(in) :: adjoint @@ -131,7 +130,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 @@ -139,7 +138,7 @@ contains k_lo = k_ln ! Fill in loss matrix - call build_loss_matrix(loss, adjoint=adjoint) + call build_loss_matrix(loss, adjoint=adjoint) ! Fill in production matrix call build_prod_matrix(prod, adjoint=adjoint) @@ -163,9 +162,8 @@ contains case(2) cmfd_linsolver => cmfd_linsolver_2g case default - message = 'Must use PETSc for more than 2 groups' - call fatal_error() - end select + cmfd_linsolver => cmfd_linsolver_ng + end select ! Set tolerances ktol = cmfd_ktol @@ -174,16 +172,19 @@ contains end subroutine init_data !=============================================================================== -! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem +! COMPUTE_ADJOINT computes a mathematical adjoint of CMFD problem !=============================================================================== subroutine compute_adjoint() use error, only: fatal_error - use global, only: message #ifdef PETSC use global, only: cmfd_write_matrices +#else + use global, only: message +#endif +#ifdef PETSC ! Transpose matrices call loss % transpose() call prod % transpose() @@ -223,7 +224,7 @@ contains iconv = .false. ! Set up tolerances - atoli = cmfd_atoli + atoli = cmfd_atoli rtoli = cmfd_rtoli toli = rtoli*100._8 @@ -279,7 +280,7 @@ contains end do - end subroutine execute_power_iter + end subroutine execute_power_iter !=============================================================================== ! WIELANDT SHIFT @@ -313,7 +314,7 @@ contains subroutine convergence(iter, innerits, iconv) - use constants, only: ONE, TINY_BIT + use constants, only: ONE, ZERO use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV @@ -328,7 +329,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)) @@ -355,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 @@ -364,7 +366,7 @@ contains integer, intent(out) :: its ! number of inner iterations integer :: g ! group index - integer :: i ! loop counter for x + integer :: i ! loop counter for x integer :: j ! loop counter for y integer :: k ! loop counter for z integer :: n ! total size of vector @@ -397,8 +399,14 @@ 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) + call tmpx % copy(x) ! Perform red/black gs iterations REDBLACK: do irb = 0,1 @@ -455,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 @@ -464,7 +473,7 @@ contains integer, intent(out) :: its ! number of inner iterations integer :: g ! group index - integer :: i ! loop counter for x + integer :: i ! loop counter for x integer :: j ! loop counter for y integer :: k ! loop counter for z integer :: n ! total size of vector @@ -509,8 +518,14 @@ 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) + call tmpx % copy(x) ! Perform red/black gs iterations REDBLACK: do irb = 0,1 @@ -587,6 +602,105 @@ 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 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 + 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 + 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 + + ! 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) + + ! 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)) + end do + do icol = didx + 1, A % get_row(irow + 1) - 1 + tmp1 = tmp1 + A % val(icol)*x % val(A % get_col(icol)) + end do + + ! Solve for new x + x1 = (b % val(irow) - tmp1)/A % val(didx) + + ! Perform overrelaxation + x % val(irow) = (ONE - w)*x % val(irow) + w*x1 + + 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 !=============================================================================== @@ -595,7 +709,7 @@ contains use global, only: cmfd, cmfd_write_matrices, current_batch - character(len=25) :: filename ! name of file to write data + character(len=25) :: filename ! name of file to write data integer :: n ! problem size ! Get problem size @@ -608,11 +722,11 @@ contains if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) end if - ! Save values + ! Save values if (adjoint_calc) then cmfd % adj_phi = phi_n % val else - cmfd % phi = phi_n % val + cmfd % phi = phi_n % val end if ! Save eigenvalue @@ -691,8 +805,8 @@ contains subroutine finalize() - ! Destroy all objects - call loss % destroy() + ! Destroy all objects + call loss % destroy() call prod % destroy() call phi_n % destroy() call phi_o % destroy() diff --git a/src/constants.F90 b/src/constants.F90 index 7d0657b695..ae9c77287a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -393,9 +393,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 0a780aeb62..b5f9c4a01e 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 @@ -329,9 +332,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' @@ -344,11 +344,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. @@ -366,7 +364,6 @@ module global ! CMFD run logicals logical :: cmfd_on = .false. - logical :: cmfd_tally_on = .true. ! CMFD display info character(len=25) :: cmfd_display = 'balance' diff --git a/src/input_xml.F90 b/src/input_xml.F90 index d21c40edc1..869f9fd63f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -265,6 +265,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) + temp_str = to_lower(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 @@ -2639,6 +2647,12 @@ contains ! Get index of mesh filter k = t % find_filter(FILTER_MESH) + ! Check to make sure mesh filter was specified + if (k == 0) then + message = "Cannot tally surface current without a mesh filter." + call fatal_error() + end if + ! Get pointer to mesh i_mesh = t % filters(k) % int_bins(1) m => meshes(i_mesh) 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 { diff --git a/src/source.F90 b/src/source.F90 index 5f928436b9..9575020e6d 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,20 @@ contains end do end if + ! Write out initial source + if (write_initial_source) then + message = 'Writing out initial source guess...' + call write_message(1) +#ifdef HDF5 + filename = trim(path_output) // 'initial_source.h5' +#else + filename = trim(path_output) // 'initial_source.binary' +#endif + call sp % file_create(filename, serial = .false.) + call sp % write_source_bank() + call sp % file_close() + end if + end subroutine initialize_source !=============================================================================== 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