crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #315 from bhermanmit/cmfd-additions

Browse source 
Extra Additions to CMFD
This commit is contained in:
Paul Romano 2014-09-17 23:42:09 -04:00
commit 47274da357
15 changed files with 332 additions and 280 deletions
Download patch file Download diff file Expand all files Collapse all files

124
docs/source/usersguide/input.rst
View file

@ -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
``<state_point>`` Element
-------------------------
@ -1303,15 +1310,6 @@ Currently, it allows users to accelerate fission source convergence during
inactive neutron batches. To run CMFD, the ``<run_cmfd>`` element in
``settings.xml`` should be set to "true".
``<active_flush>`` Element
--------------------------
The ``<active_flush>`` 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
``<begin>`` Element
-------------------
@ -1333,7 +1331,25 @@ The ``<display>`` 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
``<dhat_reset>`` Element
------------------------
The ``<dhat_reset>`` 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
``<downscatter>`` Element
-------------------------
The ``<downscatter>`` 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
``<feedback>`` Element
----------------------
@ -1344,24 +1360,16 @@ It can be turned on with "true" and off with "false".
*Default*: false
``<inactive>`` Element
----------------------
``<gauss_seidel_tolerance>`` Element
--------------------
The ``<inactive>`` 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 ``<gauss_seidel_tolerance>`` 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
``<inactive_flush>`` Element
----------------------------
The ``<inactive_flush>`` 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 ``<num_flushes>`` element.
*Defualt*: 9999
*Default*: 1.e-10 1.e-5
``<ksp_monitor>`` Element
-------------------------
@ -1370,9 +1378,16 @@ The ``<ksp_monitor>`` 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
``<ktol>`` Element
--------------------
The ``<ktol>`` element specifies the tolerance on the eigenvalue when performing
CMFD power iteration.
*Default*: 1.e-8
``<mesh>`` Element
------------------
@ -1441,14 +1456,6 @@ not impact the calculation.
*Default*: 1.0
``<num_flushes>`` Element
-------------------------
The ``<num_flushes>`` element controls the number of CMFD tally resets that
occur during inactive CMFD batches.
*Default*: 9999
``<power_monitor>`` Element
---------------------------
@ -1461,16 +1468,8 @@ This option can be turned on with "true" and turned off with "false".
-------------------------
The ``<run_adjoint>`` element can be turned on with "true" to have an adjoint
calculation be performed on the last batch when CMFD is active.
*Default*: false
``<snes_monitor>`` Element
--------------------------
The ``<snes_monitor>`` 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
``<shift>`` Element
--------------------
The ``<shfit>`` 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
``<spectral>`` Element
----------------------
The ``<spectral>`` 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
``<stol>`` Element
------------------
The ``<stol>`` element specifies the tolerance on the fission source when performing
CMFD power iteration.
*Default*: 1.e-8
``<tally_reset>`` Element
-------------------------
The ``<tally_reset>`` element contains a list of batch numbers in which CMFD tallies
should be reset.
*Default*: None
``<write_matrices>`` Element
----------------------------

159
src/cmfd_data.F90
View file

@ -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
!===============================================================================

48
src/cmfd_execute.F90
View file

@ -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
!===============================================================================

3
src/cmfd_header.F90
View file

@ -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

62
src/cmfd_input.F90
View file

@ -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)

4
src/cmfd_power_solver.F90
View file

@ -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))

136
src/cmfd_solver.F90
View file

@ -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
@ -163,8 +162,7 @@ contains
case(2)
cmfd_linsolver => cmfd_linsolver_2g
case default
message = 'Must use PETSc for more than 2 groups'
call fatal_error()
cmfd_linsolver => cmfd_linsolver_ng
end select
! Set tolerances
@ -180,10 +178,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()
@ -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
@ -397,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)
@ -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
@ -509,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)
@ -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
!===============================================================================

3
src/constants.F90
View file

@ -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

15
src/global.F90
View file

@ -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'

8
src/input_xml.F90
View file

@ -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

25
src/relaxng/cmfd.rnc
View file

@ -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 }?
}

3
src/relaxng/settings.rnc
View file

@ -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 {

15
src/source.F90
View file

@ -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
!===============================================================================

4
tests/test_cmfd_feed/cmfd.xml
View file

@ -12,7 +12,5 @@
<display> dominance </display>
<solver> power </solver>
<feedback> true </feedback>
<atoli> 1.e-15 </atoli>
<rtoli> 1.e-20 </rtoli>
<gauss_seidel_tolerance> 1.e-15 1.e-20 </gauss_seidel_tolerance>
</cmfd>

3
tests/test_cmfd_nofeed/cmfd.xml
View file

@ -12,7 +12,6 @@
<display> dominance </display>
<solver> power </solver>
<feedback> false </feedback>
<atoli> 1.e-15 </atoli>
<rtoli> 1.e-20 </rtoli>
<gauss_seidel_tolerance> 1.e-15 1.e-20 </gauss_seidel_tolerance>
</cmfd>